static const char rcsid[] = "$Id: sva.c,v 4.1 1999/11/13 22:09:00 ghiorso Exp ghiorso $"; /* MELTS Source Code: RCS $Log: sva.c,v $ MELTS Source Code: RCS Revision 4.1 1999/11/13 22:09:00 ghiorso MELTS Source Code: RCS Master Server Version (MELTS/Calc/pMELTS). MELTS Source Code: RCS MELTS Source Code: RCS Revision 4.0 1999/06/18 17:25:43 ghiorso MELTS Source Code: RCS Java MELTS v 1.1.0 Initial Check in MELTS Source Code: RCS * Revision 3.7 1997/06/21 22:49:24 ghiorso * June 1997 MELTS 3.0.x release * (prior to new entropy and regression model being introduced) * * Revision 3.6 1997/05/03 20:23:03 ghiorso * *** empty log message *** * * Revision 3.5 1997/03/27 17:03:08 ghiorso * *** empty log message *** * * Revision 3.4 1996/09/24 20:33:18 ghiorso * Version modified for OSF/1 4.0 * * Revision 3.3 1995/12/09 19:26:38 ghiorso * Interface revisions for status box and graphics display * * Revision 3.2 1995/11/01 22:40:27 ghiorso * Implementation of subsolidus options after Asimow. * Additional implementation of nepheline solid solutions. * * Revision 3.2 1995/11/01 22:40:27 ghiorso * Implementation of subsolidus options after Asimow. * Additional implementation of nepheline solid solutions. * * Revision 3.1 1995/08/18 19:15:44 ghiorso * MELTS Version 3 - Initial Entry * */ /* **++ ** FACILITY: Thread Safe Silicate Melts Crystallization Package ** ** MODULE DESCRIPTION: ** ** Singular value analysis routines (SVDCMP and SVBKSB) described in: ** ** Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T. ** (1988) Numerical Recipes in C ** Cambridge University Press, New York, 735 pages ** ** File: SVA.C **-- */ #include #include "recipes.h" #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) #define DSQR(a) ((a)*(a)) #define DMAX(a,b) ((a) > (b) ? (a) : (b)) #define IMIN(a,b) ((a) > (b) ? (b) : (a)) static double pythag(double a, double b) { 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))); } void svdcmp(double **a, int m, int n, double w[], double **v) { int flag,i,its,j,jj,k,l,nm; double anorm,c,f,g,h,s,scale,x,y,z,*rv1; rv1=dvector(1,n); g=scale=anorm=0.0; 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--) { if (i < n) { if (g) { for (j=l;j<=n;j++) 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--) { 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--) { for (its=1;its<=30;its++) { flag=1; for (l=k;l>=1;l--) { nm=l-1; if ((double)(fabs(rv1[l])+anorm) == anorm) { flag=0; break; } if ((double)(fabs(w[nm])+anorm) == anorm) break; } if (flag) { c=0.0; s=1.0; for (i=l;i<=k;i++) { f=s*rv1[i]; rv1[i]=c*rv1[i]; if ((double)(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) { if (z < 0.0) { w[k] = -z; for (j=1;j<=n;j++) v[j][k] = -v[j][k]; } break; } if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations"); x=w[l]; 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; 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; 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; } } free_dvector(rv1,1,n); } void svbksb(double **u, double w[], double **v, int m, int n, double b[], int e[], double x[]) { int jj,j,i; double s,*tmp; tmp=dvector(1,n); for (j=1;j<=n;j++) { s=0.0; if (w[j]) { for (i=1;i<=m;i++) s += (e[i]) ? u[i][j]*b[i] : 0.0; s /= w[j]; } tmp[j]=s; } for (j=1;j<=n;j++) { s=0.0; for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj]; x[j]=s; } free_dvector(tmp,1,n); } #undef NRANSI #undef SIGN #undef DMAX #undef IMIN #undef DSQR /* end-of-file SVA.C */