SVD.C

Go to the documentation of this file.

#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include <math.h>
#include <io.h>
#include "mgui.h"
#include "DAS_Spatram.h"
#include "dcl.h"
#include "dil.h"
#include "bil.h"
#include "DOAS.h"
#include "nrutil.h"
#include "Marq.h"



#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}

#define TOL 1.0e-5 //Default value for single precision and variables scaled to order unity.

void covsrt(float **covar, int ma, int ia[], int mfit)
//Expand in storage the covariance matrix covar so as to take into account parameters that are
//being held .xed.(For the latter,return zero covariances.)
{
    int i,j,k;
    float swap;
    for (i=mfit+1;i<=ma;i++)
    for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;
    k=mfit;
    for (j=ma;j>=1;j--) 
    {
        if (ia[j])
        {
            for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
            for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
            k--;
        }
    }
}





void svdfit(float x[], float y[], float sig[], int ndata, float a[], int ma,
float **u, float **v, float w[], float *chisq,
void (*funcs)(float, float [], int))
/*
Given a set of data points x[1..ndata] y[1..ndata] with individual standard deviations
sig[1..ndata],use �� 2 minimization to determine the coe .cients a[1..ma] of the .t-
ting function y = i ai �~ afunci(x )Here we solve the .tting equations using singular
value decomposition of the ndata by ma matrix as in �� 2.6.Arrays u[1..ndata][1..ma]
v[1..ma][1..ma],and w[1..ma] provide workspace on input;on output they de .ne the
singular value decomposition and can be used to obtain the covariance matrix.The pro-
gram returns values for the ma .t parameters and �� 2 chisq The user supplies a routine
funcs(x,afunc,ma) that returns the ma basis functions evaluated at x = x in the array
afunc[1..ma].
*/
{
    void svbksb(long double **u, long double w[], long double **v, int m, int n, long double b[],
        double x[]);
    void svdcmp(long double **a, int m, int n, long double w[], long double **v);
    int j,i;
    long double wmax,tmp,thresh,sum,*b,*afunc;

    b=dvector(1,ndata);
    afunc=dvector(1,ma);
    for (i=1;i<=ndata;i++) 
    { //Accumulate coefficients of the fitting matrix.
        (*funcs)(x[i],(float *)afunc,ma);
        tmp=1.0/sig[i];
        for (j=1;j<=ma;j++) u[i][j]=afunc[j]*tmp;
            b[i]=y[i]*tmp;
    }
    svdcmp(u,ndata,ma,w,v); //Singular value decomposition.
    wmax=0.0; //Edit the singular values,given TOL from the
              //#define statement,between here ...
    for (j=1;j<=ma;j++)
        if (w[j] > wmax) wmax=w[j];
            thresh=TOL*wmax;
    for (j=1;j<=ma;j++)
        if (w[j] < thresh) w[j]=0.0; //...and here.
    svbksb(u,w,v,ndata,ma,b,a);
    *chisq=0.0; //Evaluate chi-square.
    for (i=1;i<=ndata;i++)
    {
        (*funcs)(x[i],afunc,ma);
        for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];
        *chisq += (float)(tmp=(y[i]-sum)/sig[i],tmp*tmp);
    }
    free_vector(afunc,1,ma);
    free_vector(b,1,ndata);
}



void svbksb(long double **u, long double w[], long double **v,
            int m, int n, long double b[], long double x[])
/*
Solves A � X = B for a vector X where A is specified by the arrays u[1..m][1..n],w[1..n],
v[1..n][1..n] as returned by svdcmp m and n are the dimensions of a and will be equal for
square matrices.b[1..m] is the input right-hand side.x[1..n] is the output solution vector.
No input quantities are destroyed,so the routine may be called sequentially with different bs.
*/
{
    int jj,j,i;
    long double s,*tmp;
    tmp=dvector(1,n);
    for (j=1;j<=n;j++)
    { //Calculate U T B
        s=0.0;
        if (w[j])
        { //Nonzero result only if wj is nonzero.
            for (i=1;i<=m;i++) s += u[i][j]*b[i];
            s /= w[j]; //This is the divide by wj .
        }
        tmp[j]=s;
    }
    for (j=1;j<=n;j++) 
    { //Matrix multiply by V to get answer.
        s=0.0;
        for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
        x[j]=s;
    }
    free_dvector(tmp,1,n);
}



void svdcmp(long double **a, int m, int n, long double w[], long double **v)
/*
Given a matrix a[1..m][1..n],this routine computes its singular value decomposition,
A =U � W � V(Transposed).
The matrix U replaces a on output.The diagonal matrix of singular values W
 is output as a vector w[1..n].The matrix V (not the transpose V T )is output as v[1..n][1..n].
*/
{
    FILE *fpp;
    long double pythag(long double a, long double b);
    int flag,i,its,j,jj,k,l,nm;
    long double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
    rv1=dvector(1,n);
    g=scale=anorm=0.0; //Householder reduction to bidiagonal form.
    for (i=1;i<=n;i++)
    {
        l=i+1;
        rv1[i]=scale*g;
        g=s=scale=0.0;
        if (i <= m)
        {
            for (k=i;k<=m;k++)
                scale += fabs(a[k][i]);
            if (scale)
            {
                for (k=i;k<=m;k++)
                {
                    a[k][i] /= scale;
                    s += a[k][i]*a[k][i];
                }
                f=a[i][i];
                g = -SIGN(sqrt(s),f);
                h=f*g-s;
                a[i][i]=f-g;
                for (j=l;j<=n;j++)
                {
                    for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
                        f=s/h;
                    for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
                }
                for (k=i;k<=m;k++) a[k][i] *= scale;
            }
        }
        w[i]=scale *g;
        g=s=scale=0.0;
        if (i <= m && i != n)
        {
            for (k=l;k<=n;k++) scale += fabs(a[i][k]);
            if (scale)
            {
                for (k=l;k<=n;k++) 
                {
                    a[i][k] /= scale;
                    s += a[i][k]*a[i][k];
                }
                f=a[i][l];
                g = -SIGN(sqrt(s),f);
                h=f*g-s;
                a[i][l]=f-g;
                for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
                for (j=l;j<=m;j++) 
                {
                    for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
                    for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
                }
                for (k=l;k<=n;k++) a[i][k] *= scale;
            }
        }
        anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
    }
    for (i=n;i>=1;i--)
    { //Accumulation of right-hand transformations.
        if (i < n)
        {
            if (g)
            {
                for (j=l;j<=n;j++) //Double division to avoid possible under .o .
                    v[j][i]=(a[i][j]/a[i][l])/g;
                for (j=l;j<=n;j++)
                {
                    for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
                    for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
                }
            }
            for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
        }
        v[i][i]=1.0;
        g=rv1[i];
        l=i;
    }
    for (i=IMIN(m,n);i>=1;i--) 
    { //Accumulation of left-hand transformations.
        l=i+1;
        g=w[i];
        for (j=l;j<=n;j++) a[i][j]=0.0;
        if (g)
        {
            g=1.0/g;
            for (j=l;j<=n;j++)
            {
                for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
                    f=(s/a[i][i])*g;
                for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
            }
            for (j=i;j<=m;j++) a[j][i] *= g;
        }
        else 
            for (j=i;j<=m;j++) a[j][i]=0.0;
                ++a[i][i];
    }

    for (k=n;k>=1;k--) 
    { //Diagonalization of the bidiagonal form:Loop over singular values,and over allo ed iterations.
        for (its=1;its<=30;its++) 
        {
            flag=1;
            for (l=k;l>=1;l--) 
            { //Test for splitting.
                nm=l-1; //Note that rv1[1] is always zero.
                if ((float)(fabs(rv1[l])+anorm) == anorm) 
                {
                    flag=0;
                    break;
                }
                if ((float)(fabs(w[nm])+anorm) == anorm) break;
            }
            if (flag)
            {
                c=0.0; //Cancellation of rv1[l],if l > 1
                s=1.0;
                for (i=l;i<=k;i++)
                {
                    f=s*rv1[i];
                    rv1[i]=c*rv1[i];
                    if ((float)(fabs(f)+anorm) == anorm) break;
                    g=w[i];
                    h=pythag(f,g);
                    w[i]=h;
                    h=1.0/h;
                    c=g*h;
                    s = -f*h;
                    for (j=1;j<=m;j++)
                    {
                        y=a[j][nm];
                        z=a[j][i];
                        a[j][nm]=y*c+z*s;
                        a[j][i]=z*c-y*s;
                    }
                }
            }
            z=w[k];
            if (l == k)
            { //Convergence.
                if (z < 0.0)
                { //Singular value is made nonnegative.
                    w[k] = -z;
                    for (j=1;j<=n;j++) v[j][k] = -v[j][k];
                }
                break;
            }
            if (its == 30) nrerror("no convergence in 30 svdcmp iterations");
            x=w[l]; //Shift from bottom 2-by-2 minor.
            nm=k-1;
            y=w[nm];
            g=rv1[nm];
            h=rv1[k];
            f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
            g=pythag(f,1.0);
            f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
            c=s=1.0; //Next QR transformation:
            for (j=l;j<=nm;j++)
            {
                i=j+1;
                g=rv1[i];
                y=w[i];
                h=s*g;
                g=c*g;
                z=pythag(f,h);
                rv1[j]=z;
                c=f/z;
                s=h/z;
                f=x*c+g*s;
                g = g*c-x*s;
                h=y*s;
                y *=c;
                for (jj=1;jj<=n;jj++)
                {
                    x=v[jj][j];
                    z=v[jj][i];
                    v[jj][j]=x*c+z*s;
                    v[jj][i]=z*c-x*s;
                }
                z=pythag(f,h);
                w[j]=z; //Rotation can be arbitrary if z =0
                if (z)
                {
                    z=1.0/z;
                    c=f*z;
                    s=h*z;
                }
                f=c*g+s*y;
                x=c*y-s*g;
                for (jj=1;jj<=m;jj++)
                {
                    y=a[jj][j];
                    z=a[jj][i];
                    a[jj][j]=y*c+z*s;
                    a[jj][i]=z*c-y*s;
                }
            }
            rv1[l]=0.0;
            rv1[k]=f;
            w[k]=x;
        }
    }

            fpp = fopen("matUInSVD.dat", "w");
        for(k=1;k<=m;k++)
        {
            for(l=1;l<=m;l++)
            {
                fprintf(fpp,"%7.6e ",a[k][l]) ;
            }
            fprintf(fpp,"\n");
        }

        fclose(fpp);


    free_vector(rv1,1,n);
}



long double pythag(long double a, long double b)
//Computes (a 2 + b 2 ) 1/2 without destructive underflow or overflow.
{
    long double absa,absb;
    absa=fabs(a);
    absb=fabs(b);
    if (absa > absb) return absa*sqrt(1.0+DSQR(absb/absa));
    else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb)));
}