static const char rcsid[] = "$Id: nrutil.c,v 4.1 1999/11/13 22:09:00 ghiorso Exp ghiorso $"; /* MELTS Source Code: RCS $Log: nrutil.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.6 1997/06/21 22:49:36 ghiorso * June 1997 MELTS 3.0.x release * (prior to new entropy and regression model being introduced) * * Revision 3.5 1997/05/03 20:23:12 ghiorso * *** empty log message *** * * Revision 3.4 1997/03/27 17:03:20 ghiorso * *** empty log message *** * * Revision 3.3 1996/09/24 20:33:27 ghiorso * Version modified for OSF/1 4.0 * * Revision 3.2 1995/12/09 19:26:38 ghiorso * Interface revisions for status box and graphics display * * Revision 3.1 1995/08/18 18:03:34 ghiorso * MELTS Version 3 - Initial Entry * * Revision 3.1 1995/08/18 18:03:34 ghiorso * MELTS Version 3 - Initial Entry * */ /* **++ ** FACILITY: Thread Safe Silicate Melts Crystallization Package ** ** MODULE DESCRIPTION: ** ** Support routines 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: NRUTIL.C **-- */ #include #include #include #include "recipes.h" void nrerror(error_text) char error_text[]; { void abort(); fprintf(stderr,"Numerical Recipes run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); abort(); } double *vector(nl,nh) int nl,nh; { double *v; v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double)); if (!v) nrerror("allocation failure in vector()"); return v-nl; } int *ivector(nl,nh) int nl,nh; { int *v; v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int)); if (!v) nrerror("allocation failure in ivector()"); return v-nl; } double *dvector(nl,nh) int nl,nh; { double *v; v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double)); if (!v) nrerror("allocation failure in dvector()"); return v-nl; } double **matrix(nrl,nrh,ncl,nch) int nrl,nrh,ncl,nch; { int i; double **m; m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*)); if (!m) nrerror("allocation failure 1 in matrix()"); m -= nrl; for(i=nrl;i<=nrh;i++) { m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double)); if (!m[i]) nrerror("allocation failure 2 in matrix()"); m[i] -= ncl; } return m; } double **dmatrix(nrl,nrh,ncl,nch) int nrl,nrh,ncl,nch; { int i; double **m; m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*)); if (!m) nrerror("allocation failure 1 in dmatrix()"); m -= nrl; for(i=nrl;i<=nrh;i++) { m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double)); if (!m[i]) nrerror("allocation failure 2 in dmatrix()"); m[i] -= ncl; } return m; } int **imatrix(nrl,nrh,ncl,nch) int nrl,nrh,ncl,nch; { int i,**m; m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*)); if (!m) nrerror("allocation failure 1 in imatrix()"); m -= nrl; for(i=nrl;i<=nrh;i++) { m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int)); if (!m[i]) nrerror("allocation failure 2 in imatrix()"); m[i] -= ncl; } return m; } double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl) double **a; int oldrl,oldrh,oldcl,oldch,newrl,newcl; { int i,j; double **m; m=(double **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(double*)); if (!m) nrerror("allocation failure in submatrix()"); m -= newrl; for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl; return m; } void free_vector(v,nl,nh) double *v; int nl,nh; { free((void*) (v+nl)); } void free_ivector(v,nl,nh) int *v,nl,nh; { free((void*) (v+nl)); } void free_dvector(v,nl,nh) double *v; int nl,nh; { free((void*) (v+nl)); } void free_matrix(m,nrl,nrh,ncl,nch) double **m; int nrl,nrh,ncl,nch; { int i; for(i=nrh;i>=nrl;i--) free((void*) (m[i]+ncl)); free((void*) (m+nrl)); } void free_dmatrix(m,nrl,nrh,ncl,nch) double **m; int nrl,nrh,ncl,nch; { int i; for(i=nrh;i>=nrl;i--) free((void*) (m[i]+ncl)); free((void*) (m+nrl)); } void free_imatrix(m,nrl,nrh,ncl,nch) int **m; int nrl,nrh,ncl,nch; { int i; for(i=nrh;i>=nrl;i--) free((void*) (m[i]+ncl)); free((void*) (m+nrl)); } void free_submatrix(b,nrl,nrh,ncl,nch) double **b; int nrl,nrh,ncl,nch; { free((void*) (b+nrl)); } double **convert_matrix(a,nrl,nrh,ncl,nch) double *a; int nrl,nrh,ncl,nch; { int i,j,nrow,ncol; double **m; nrow=nrh-nrl+1; ncol=nch-ncl+1; m = (double **) malloc((unsigned) (nrow)*sizeof(double*)); if (!m) nrerror("allocation failure in convert_matrix()"); m -= nrl; for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl; return m; } void free_convert_matrix(b,nrl,nrh,ncl,nch) double **b; int nrl,nrh,ncl,nch; { free((void*) (b+nrl)); } #define ITMAX 100 #define EPS 3.0e-10 double zbrent(double (*func)(double),double xx1,double xx2,double tol) { int iter; double a=xx1,b=xx2,c,d,e,min1,min2; double fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm; void nrerror(); if (fb*fa > 0.0) { printf("Warning: Root not bracketed in ZBRENT (%g,%g)\n", fa, fb); if (fabs(fa) < fabs(fb)) return a; else return b; } fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (fabs(fc) < fabs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*fabs(b)+0.5*tol; xm=0.5*(c-b); if (fabs(xm) <= tol1 || fb == 0.0) return b; if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=fabs(p); min1=3.0*xm*q-fabs(tol1*q); min2=fabs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (fabs(d) > tol1) b += d; else b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1)); fb=(*func)(b); } nrerror("Maximum number of iterations exceeded in ZBRENT"); } #undef ITMAX #undef EPS /* end of file NRUTIL.C */