Gram-Schmidt Method

gram-schmidt.c

/**
 * Written by Andy Pan
 * Language: C
 * Algorithm: Gram-Schmidt Method
 * Creation Date: 2014/12/6
 * Last Modified: 2014/12/7
**/
#include<stdio.h>
#include<math.h>

void inputMatrix(double* matrix,int m,int n){
  int i,c=m*n;
  for(i=0;i<c;i++)
    scanf("%lf",matrix+i);
}

void transposeMatrix(double* matrix,int m,int n,double* mat2){
  int i,j;
  for(i=0;i<n;i++)
    for(j=0;j<m;j++)
      mat2[i*m+j]=matrix[j*n+i];
}

void dumpMatrix(double* matrix,int m,int n){
  int i,j,c=0;
  for(i=0;i<m;i++){
    for(j=0;j<n;j++){
      //this format is correspondent to the example program for comparing
      printf("%10.4lf ",matrix[c]);
      c++;
    }
    printf("\n");
  }
}

void multiplyMatrix(double* a,double* b,int m,int n,int p,double* r){
  int i,j,k;
  for(i=0;i<m;i++){
    for(j=0;j<n;j++){
      double v=0;
      for(k=0;k<n;k++)
        v+=a[i*n+k]*b[k*p+j];
      r[i*p+j]=v;
    }
  }
}

double dot(double a[],double b[],int dimen){
  int i;
  double ans=0;
  for(i=0;i<dimen;i++)
    ans+=a[i]*b[i];
  return ans;
}

double length2(double a[],int dimen){
  //length square=a dot a
  return dot(a,a,dimen);
}

void add(double a[],double b[],int dimen,double r[]){
  int i;
  for(i=0;i<dimen;i++)
    r[i]=a[i]+b[i];
}

void multiply(double x[],double a,int dimen,double r[]){
  int i;
  for(i=0;i<dimen;i++)
    r[i]=x[i]*a;
}

void clone(double a[],int dimen,double r[]){
  int i;
  for(i=0;i<dimen;i++)
    r[i]=a[i];
}

void proj(double a[],double e[],int dimen,double r[]){
  double q=dot(e,a,dimen)/dot(e,e,dimen);
  multiply(e,q,dimen,r);
}

void normalize(double x[],int dimen,double r[]){
  double len=sqrt(length2(x,dimen));
  multiply(x,1.0/len,dimen,r);
  //cf. round-off error: x*(1.0/a) vs. x/a
  //int i;
  //for(i=0;i<dimen;i++)
  //  r[i]=x[i]/len;
}

void QRComputeQ(double* x,int m,int n,double* Q){
  int i,j;
  double *ai,*ui,*ej;
  double u[m][n],tmp[n];
  for(i=0;i<m;i++){
    ai=x+i*n;
    ui=u[i];
    //u_m = a_m - sigma [i=1 to m-1] ( proj_(e_i) a_m)
    clone(ai,n,ui);
    for(j=0;j<i;j++){
      ej=Q+j*n;
      proj    (ai,ej,n,tmp);
      multiply(tmp,-1,n,tmp);
      add     (ui,tmp,n,ui);
    }
    //Q=[e1,e2,e3,...,e_n] so we save it directly
    normalize(ui,n,Q+i*n);
  }
}

void QRComputeR(double* x,double* Q,int m,int n,double* R){
  int i,j,c=0;
  for(i=0;i<m;i++){
    for(j=0;j<m;j++){
      if(j>=i){
        R[c]=dot(Q+i*n,x+j*n,n);
      }else{
        R[c]=0;
      }
      c++;
    }
  }
}

void QRDecomposition(double* A,int m,int n,double* Q,double* R){
  double AT[n][m],QT[n][m];
  transposeMatrix( A     ,m,n,*AT);
  QRComputeQ     (*AT    ,n,m,*QT);
  transposeMatrix(*QT    ,n,m, Q );
  QRComputeR     (*AT,*QT,n,m, R );
  
}

double QInftyNorm(double* Q,int m,int n){
  int i,j;
  double rtn=0,tmp,sum,QT[n][m],QTQ[n][n];
  transposeMatrix(Q,m,n,*QT);
  multiplyMatrix(*QT,Q,n,m,n,*QTQ);
  //infinity norm is defined as the max absolute of row sum
  for(i=0;i<n;i++){
    sum=0;
    //(Q^T)Q-I, ad-hoc
    for(j=0;j<n;j++){
      tmp=QTQ[i][j]-(i==j?1:0);
      sum+=tmp>0?tmp:-tmp;
    }
    rtn=(sum>rtn)?sum:rtn;
  }
  return rtn;
}

void printMartixComplete(double* x,int m,int n,const char* STR){
  printf(STR);
  printf("[%d,%d]=\n",m,n);
  dumpMatrix(x,m,n);
  //optional: add an extra line at the end
  //printf("\n");
}

int main(){
  int m,n;
  while(scanf("%d%d",&m,&n)){
  
    if(m==0&&n==0) return 0;
    
    double matA   [m][n],
         matQ1  [m][n],
         matR1  [n][n],
         matQ1R1[m][n];
         
    inputMatrix(*matA,m,n);
    
    QRDecomposition(*matA,m,n,*matQ1,*matR1);
    
    multiplyMatrix(*matQ1,*matR1,m,n,n,*matQ1R1);
    
    printMartixComplete(*matA,m,n,"A");
    printMartixComplete(*matQ1,m,n,"Q1");
    printf("Orthogonality of Q1 = %.14le\n",QInftyNorm(*matQ1,m,n));
    
    printMartixComplete(*matR1,n,n,"R1");
    printMartixComplete(*matQ1R1,m,n,"Q1R1");
    
    printf("===========================================\n");
  }
  return 0;
}