matrixutil.cpp

Go to the documentation of this file.

#include "matrixutil.hpp"
#include <cmath>
#include <iostream>
#include "vc6compat.h"
using namespace std;

#define TRACE(a) 
#define TRACEX(a, b)
#define SQR(X) ((X)*(X))


bool QR (int zeilen, int spalten, double **A, double *d, double *b, const double epsilon) 
{ 
 if (zeilen<spalten) {
  throw("Matrix darf nicht weniger Zeilen als Spalten haben !");
 } 
 //
 bool singular=false;
 double altenorm,hilfsnorm,beta, produkt;
 altenorm=0;
 for (int j=0; j<spalten; j++) {
  // jede Spalte muss transformiert werden
  // hier muss als erstes ||a|| berechnet werden
  // dazu brauche ich eine "Hilfsnorm"
  hilfsnorm=0;
  for (int my=j+1; my<zeilen; my++) {
    hilfsnorm+=SQR(A[my][j]);
  }
  // jetzt ist sqrt(hilfsnorm+SQR(A[j][j]) die 2-Norm von a
  // berechne den Vektor b dorthin, wo a war
  // wenn das erste Element viel kleiner ist als die Hilfsnorm
  // dann kann man es als Null ansehen
  TRACEX("Vor der Transformation :",A);
  TRACE(hilfsnorm);
  
  if (SQR(A[j][j])<SQR(epsilon)*hilfsnorm) {
     A[j][j]=sqrt(hilfsnorm);
     d[j]=-A[j][j];
     cerr<<"Diagonal-elem. "<<j<<" = 0"<<endl;
  } else {
     d[j]=-signum(A[j][j])*sqrt(hilfsnorm+SQR(A[j][j]));
     // speicher k in d_j

     A[j][j]-=d[j];
  }
  // ein naives Abbruchkriterium:
  // Wenn das aktuelle k sehr viel kleiner ist als das letzte,
  // dann sieh die Matrix als singulär an
  // ich hab sonst nix gescheites gefunden
  // auch MATLAB und octave berechnen anstandslos die
  // QR-Zerlegung von singulären Matrizen
  
  if (fabs(d[j])<epsilon*altenorm) {
    singular=true;
    break;
  } 

  altenorm=fabs(d[j]);
  
  // berechne beta
  beta=2.0/(hilfsnorm+SQR(A[j][j]));
  
  TRACEX("nach der Transformation "<<j,A);
  TRACE(beta); 
  // jetzt sind alle Zutaten beieinander, wir können anfangen
  // den Rest der Matrix zu transformieren
  
  for (int k=j+1; k<spalten; k++) {
    // jede Spalte weiter rechts
    // berechne das innere Produkt aus b und der Spalte ak
    produkt=0;
    for (int z=j; z<zeilen; z++) {
      produkt+=A[z][j]*A[z][k];
    }
    produkt*=beta;
    TRACE(produkt);
    // hetze die Transformation auf jede Spalte a_k
    for (int z=j; z<zeilen; z++) {
      A[z][k]-=produkt*A[z][j];
    }  
  }
  // jetzt das gleiche nochmal mit der rechten Seite
  produkt=0;
  for (int z=j; z<zeilen; z++) {
    produkt+=A[z][j]*b[z];
  }
  produkt*=beta;
  TRACE(produkt);
  // hetze die Transformation auf die rechte Seite
  for (int z=j; z<zeilen; z++) {
    b[z]-=produkt*A[z][j];
  }
 
  TRACE(A);
 }
 if (singular) return true;
 return false;
} 


void rdsolve(int spalten, double **R, double *d, double *c, double *x)
{
  double tempx;
  for (int j=spalten-1; j>=0; j--) {
   tempx=c[j];
   // umgeformte rechte Seite in tempx speichern
   // abziehen die bekannten Variablen
   for (int k=j+1; k<spalten; k++)
     tempx-=R[j][k]*x[k];

   x[j]=tempx/d[j];
  } // for j
} // rsolve



// Berechne Kovarianzmatrix, aus dem Source von Gnuplot
// by Carsten Grammers, public domain
//
void Invert_RtR(double **R, double **I, int n)
{
    int i, j, k;

    /* fill in the I matrix, and check R for regularity : */

    for (i = 0; i < n; i++) {
        for (j = 0; j < i; j++) /* upper triangle isn't needed */
            I[i][j] = 0;
        I[i][i] = 1;
        if (!R[i][i])
            throw("Singular matrix in Invert_RtR");
    }

    /* Forward substitution: Solve R^T * B = I, store B in place of I */

    for (k = 0; k < n; k++) {
        for (i = k; i < n; i++) {       /* upper half needn't be computed */
            double s = I[i][k];
            for (j = k; j < i; j++)     /* for j<k, I[j][k] always stays zero! */
                s -= R[j][i] * I[j][k];
            I[i][k] = s / R[i][i];
        }
    }
    /* Backward substitution: Solve R * A = B, store A in place of B */

    for (k = 0; k < n; k++) {
        for (i = n - 1; i >= k; i--) {  /* don't compute upper triangle of A */
            double s = I[i][k];
            for (j = i + 1; j < n; j++)
                s -= R[i][j] * I[j][k];
            I[i][k] = s / R[i][i];
        }
    }
}

void compute_errors(int num_param, int num_data, double **R, double *d, double *b, double *errors) 
{
  // berechne reduziertes chi^2
  // Der Fehler steckt im Residualvektor der QR-Zerlegung
  double chisq=0;
  for (int i=num_param; i<num_data; i++) {
     chisq+=SQR(b[i]);
  }

  // geht nur für num_param<num_data - klar
  if (num_param>=num_data) throw("Error: Parameter errors only für num_data < num_data");
  chisq=sqrt(chisq/(num_data-num_param)); 
  
   // Setze Diagonalelemente in R ein
  for (int r=0; r<num_param; r++) {
    R[r][r]=d[r];
  }
 
  double **C;
  reserve_2D(C, num_param, num_param);
  Invert_RtR(R, C, num_param);
  
  for (int r=0; r<num_param; r++) {
    errors[r]=sqrt(C[r][r])*chisq;
  }
 
  delete_2D(C);
  // kovarianzmatrix löschen
}  

---