inversa1.cpp - #include#include#include#include\"nrutil.h <math.h> <stdlib.h> <stdio.h>#define LENGTH 400#define TINY 1.0e-20#define NULL1 0 void

# inversa1.cpp - #include#include#include#include"nrutil.h...

• 3

This preview shows page 1 - 3 out of 3 pages.

#include "nrutil.h" #include <math.h> #include <stdlib.h> #include <stdio.h> #define LENGTH 400 #define TINY 1.0e-20 #define NULL1 0 void ludcmp1(float **m, int n, int *indx, float *d) { int i, imax, j, kkk; float big, dum, sum, temp; float *vv, *v1; vv= new float [n]; v1= new float [n]; if(vv == NULL1 || v1 == NULL1) { printf("Can't allocation memory in ludcmp1 !"); exit(1); } *d = 1.0; for(i=0; i<n; i++){ big=0.0; for(j=0;j<n;j++) if ((temp=fabs(m[i][j]))>big) big=temp; if(big == 0.0) return; vv[i]=1.0/big; } for(j=0; j<n; j++){ for(i=0; i<j; i++){ sum=m[i][j]; for(kkk=0; kkk<i; kkk++) sum -= m[i][kkk]*m[kkk][j]; m[i][j]=sum; } big=0.0; for (i=j; i<n; i++){ sum=m[i][j]; for (kkk=0; kkk<j; kkk++) sum -= m[i][kkk]*m[kkk][j]; m[i][j]=sum; if ((dum=vv[i]*fabs(sum)) >= big){ big=dum; imax=i; } } if (j != imax){ for(kkk=0; kkk<n; kkk++){ dum=m[imax][kkk]; m[imax][kkk]=m[j][kkk]; m[j][kkk]=dum; } *d = -(*d); vv[imax]=vv[j]; } indx[j]=imax; if (fabs(m[j][j])<TINY){ printf("divided by zero!\n");
m[j][j]=TINY; } if (j != n){ dum=1.0/(m[j][j]); for(i=j+1; i<n; i++) m[i][j] *= dum; } } delete [] vv; vv = NULL1; delete [] v1; v1 = NULL1; } void lubksb1(float **m, int n, int *indx, float *b) { int i, ii=-1, ip, j; float sum; for (i=0; i<n; i++){ ip=indx[i]; sum=b[ip]; b[ip]=b[i]; if (ii>=0) for (j=ii; j<i; j++) sum -= m[i][j]*b[j];