00001 #include "vc6compat.h"
00002 #include "matrixutil.hpp"
00003 #include <cmath>
00004 using namespace std;
00005
00006 #define TRACE(a)
00007 #define TRACEX(a, b)
00008 #define SQR(X) ((X)*(X))
00009
00010
00011 bool QR (int zeilen, int spalten, double **A, double *d, double *b, const double epsilon)
00012 {
00013 if (zeilen<spalten) {
00014 throw("Matrix darf nicht weniger Zeilen als Spalten haben !");
00015 }
00016
00017 bool singular=false;
00018 double altenorm,hilfsnorm,beta, produkt;
00019 altenorm=0;
00020 for (int j=0; j<spalten; j++) {
00021
00022
00023
00024 hilfsnorm=0;
00025 for (int my=j+1; my<zeilen; my++) {
00026 hilfsnorm+=SQR(A[my][j]);
00027 }
00028
00029
00030
00031
00032 TRACEX("Vor der Transformation :",A);
00033 TRACE(hilfsnorm);
00034
00035 if (SQR(A[j][j])<SQR(epsilon)*hilfsnorm) {
00036 A[j][j]=sqrt(hilfsnorm);
00037 d[j]=A[j][j];
00038 } else {
00039 d[j]=-signum(A[j][j])*sqrt(hilfsnorm+SQR(A[j][j]));
00040
00041
00042 A[j][j]-=d[j];
00043 }
00044
00045
00046
00047
00048
00049
00050
00051 if (fabs(d[j])<epsilon*altenorm) {
00052 singular=true;
00053 break;
00054 }
00055
00056 altenorm=fabs(d[j]);
00057
00058
00059 beta=2.0/(hilfsnorm+SQR(A[j][j]));
00060
00061 TRACEX("nach der Transformation "<<j,A);
00062 TRACE(beta);
00063
00064
00065
00066 for (int k=j+1; k<spalten; k++) {
00067
00068
00069 produkt=0;
00070 for (int z=j; z<zeilen; z++) {
00071 produkt+=A[z][j]*A[z][k];
00072 }
00073 produkt*=beta;
00074 TRACE(produkt);
00075
00076 for (int z=j; z<zeilen; z++) {
00077 A[z][k]-=produkt*A[z][j];
00078 }
00079 }
00080
00081 produkt=0;
00082 for (int z=j; z<zeilen; z++) {
00083 produkt+=A[z][j]*b[z];
00084 }
00085 produkt*=beta;
00086 TRACE(produkt);
00087
00088 for (int z=j; z<zeilen; z++) {
00089 b[z]-=produkt*A[z][j];
00090 }
00091
00092 TRACE(A);
00093 }
00094 if (singular) return true;
00095 return false;
00096 }
00097
00104
00105 void rdsolve(int spalten, double **R, double *d, double *c, double *x)
00106 {
00107 double tempx;
00108 for (int j=spalten-1; j>=0; j--) {
00109 tempx=c[j];
00110
00111
00112 for (int k=j+1; k<spalten; k++)
00113 tempx-=R[j][k]*x[k];
00114
00115 x[j]=tempx/d[j];
00116 }
00117 }
00118
00119
00120
00121
00122
00123
00124 void Invert_RtR(double **R, double **I, int n)
00125 {
00126 int i, j, k;
00127
00128
00129
00130 for (i = 0; i < n; i++) {
00131 for (j = 0; j < i; j++)
00132 I[i][j] = 0;
00133 I[i][i] = 1;
00134 if (!R[i][i])
00135 throw("Singular matrix in Invert_RtR");
00136 }
00137
00138
00139
00140 for (k = 0; k < n; k++) {
00141 for (i = k; i < n; i++) {
00142 double s = I[i][k];
00143 for (j = k; j < i; j++)
00144 s -= R[j][i] * I[j][k];
00145 I[i][k] = s / R[i][i];
00146 }
00147 }
00148
00149
00150 for (k = 0; k < n; k++) {
00151 for (i = n - 1; i >= k; i--) {
00152 double s = I[i][k];
00153 for (j = i + 1; j < n; j++)
00154 s -= R[i][j] * I[j][k];
00155 I[i][k] = s / R[i][i];
00156 }
00157 }
00158 }
00159
00160 void compute_errors(int num_param, int num_data, double **R, double *d, double *b, double *errors)
00161 {
00162
00163
00164 double chisq=0;
00165 for (int i=num_param; i<num_data; i++) {
00166 chisq+=SQR(b[i]);
00167 }
00168
00169
00170 if (num_param>=num_data) throw("Error: Parameter errors only für num_data < num_data");
00171 chisq=sqrt(chisq/(num_data-num_param));
00172
00173
00174 for (int r=0; r<num_param; r++) {
00175 R[r][r]=d[r];
00176 }
00177
00178 double **C;
00179 reserve_2D(C, num_param, num_param);
00180 Invert_RtR(R, C, num_param);
00181
00182 for (int r=0; r<num_param; r++) {
00183 errors[r]=sqrt(C[r][r])*chisq;
00184 }
00185
00186 delete_2D(C);
00187
00188 }