Category : Science and Education
Archive   : LAPLACE.ZIP
Filename : SOLVE.C

 
Output of file : SOLVE.C contained in archive : LAPLACE.ZIP
/* decomp & solve - solution of linear system


AUTHORS
Forsythe, Malcolm, and Moler, from "Computer Methods for
Mathematical Computations"

translated to C by J. R. Van Zandt

#define TEST
*/

#include
#include

#define aa(i,j) a[ix[i]+j] /* like a[i][j] or a[i*ndim+j] but faster */

double decomp(ndim,n,a,ipvt) /* returns estimate of condition number */
int ndim, /* 2nd declared dimension of a (# columns) */
n, /* rows 0...n-1 and columns 0...n-1 of a are used */
ipvt[]; /* pivot columns */
double a[]; /* two dimensional array of matrix coefficients */
{ double ek, t, anorm, ynorm, znorm, cond;
int nm1, i, j, k, kp1, kb, km1, m;
static double work[25];
static int ix[25];

if(n>25)
{fprintf(stderr,"decomp: matrix too big\n");
exit(1);
}
/*
double *work;
int *ix;
work=malloc(n*sizeof(double));
ix=malloc(n*sizeof(int));
if(work==0 || ix==0)
{fprintf(stderr,"decomp: no workspace\n");
exit(1);
}
*/
for (i=0; i ix[i]=i*ndim;


nm1=n-1;
ipvt[nm1]=0;
if(n==1)
{
/*
free(work);
free(ix);
*/
if(a[0]==0.)
return 1.e32;
else
return 1.;
}
/* compute 1-norm of a */
anorm=0.;
for (j=0; j {t=0.;
for (i=0; i t += fabs(aa(i,j));
if(t>anorm)
anorm=t;
}
/* Gaussian elimination with partial pivoting */
for (k=0; k {kp1=k+1;
/* Find pivot */
m=k;
for (i=kp1; i if(fabs(aa(i,k)) > fabs(aa(m,k)))
m=i;
ipvt[k]=m;
t=aa(m,k);
/* printf("pivoting on a[%d][%d] = %8.4f\n",m,k,t); *****/
if(m!=k)
{ipvt[nm1]=-ipvt[nm1];
aa(m,k)=aa(k,k);
aa(k,k)=t;
}
/* skip step if pivot is zero */
if(t!=0.)
{ /* compute multipliers */
for (i=kp1; i aa(i,k) = -aa(i,k)/t;
/* interchange and eliminate by columns */
for (j=kp1; j {t=aa(m,j);
aa(m,j) = aa(k,j);
aa(k,j) = t;
if(t != 0.)
for (i=kp1; i aa(i,j) += aa(i,k)*t;
}
}
/* showm("After pivoting",a,ndim,n); /**/
}
/*
cond = (1-norm of a)*(an estimate of 1-norm of
a-inverse) estimate obtained by one step of
inverse iteration for the small singular
vector. This involves solving two systems of
equations, (a-transpose)*y = e and a*z=y there
e is a vector of +1 or -1 chosen to cause
growth in y. Estimate = (1-norm of z)/(1-norm
of y)
*/

/* solve (a-transpose)*y - e */
/* printf("5");
/* showiv(" pivot cols:", ipvt, n); /**/
for (k=0; k {t=0.;
if(k)
for (i=0; i t += aa(i,k)*work[i];
if(t<0.) ek = -1.; else ek=1.;
if (aa(k,k) == 0.)
{
/*
free(work);
free(ix);
*/
return 1.e32;
}
work[k] = -(ek+t)/aa(k,k);
}
/* printf("6");
/* showiv(" pivot cols:", ipvt, n); /**/
/* showv("decompose: work1",work,n); /**/
for (k=n-2; k>=0; k--)
{t=0.;
for (i=k+1; i t += aa(i,k)*work[k];

work[k] = t;
m=ipvt[k];
if (m!=k)
{t=work[m];
work[m] = work[k];
work[k]=t;
}
}
/* printf("8");
/* showiv(" pivot cols:", ipvt, n); /**/
ynorm=0.;
for (i=0; i ynorm += fabs(work[i]);
/* solve a*z=y */
/* printf("9");
/* showiv(" pivot cols:", ipvt, n); /**/
solve(ndim, n, a, work, ipvt);
/* printf("0"); /**/
znorm=0.;
for (i=0; i znorm += fabs(work[i]);
/* estimate condition */
cond=anorm*znorm/ynorm;
if(cond<1.)
cond=1.;
/*
free(work);
free(ix);
*/
/* printf("1"); */
return cond;
}

double invert(ndim,n,a,x) /* returns estimate of condition number of matrix */
int ndim, /* 2nd declared dimension of a and x (# columns) */
n; /* rows 0...n-1 and columns 0...n-1 of a are used */
double a[], /* matrix to be inverted */
x[]; /* resulting inverse */
{ int i, j;
double cond, condp1;
static double work[25];
static int ipvt[25];
if(n>25)
{fprintf(stderr,"invert: matrix too big\n");
exit(1);
}
/*
double *work;
int *ipvt;
ipvt=malloc(n*sizeof(int));
work=malloc(n*sizeof(double));
if(ipvt==0 || work==0)
{fprintf(stderr,"invert: not enough memory\n");
exit(1);
}
*/
cond=decomp(ndim,n,a,ipvt);
/* printf("invert: cond = %f\n\n",cond);
printf("invert: ipvt\n");
for (i=0; i printf("%8d \n",ipvt[i]);
*/
condp1=cond+1.;
if(condp1==cond)
{
/*
free(work);
free(ipvt);
*/
return cond;
}
for (i=0; i {for (j=0; j work[j]=0.;
work[i]=1.;
/* showv("invert: RHS",work,n); */
solve(ndim,n,a,work,ipvt);
for (j=0; j x[j*ndim+i]=work[j];
}
/*
free(work);
free(ipvt);
*/
return cond;
}

solve (ndim,n,a,b,ipvt)
int ndim, /* 2nd declared dimension of a (# columns) */
n, /* order of matrix a: rows 0...n-1 & columns 0...n-1 are used */
ipvt[]; /* pivot vector obtained from decomp */
double a[], /* triangularized matrix obtained from decomp */
b[]; /* right hand side vector */
{ int kb,km1,nm1,kp1,i,k,m;
double t;
static int ix[25];
if(n>25)
{fprintf(stderr,"solve: matrix too big\n");
exit(1);
}
/*
int *ix;
ix=malloc(n*sizeof(int));
if(ix==0)
{fprintf(stderr,"solve: no workspace\n");
exit(1);
}
*/
/* printf("a"); */
for (i=0; i ix[i]=i*ndim;
/* showm("solve: decomposed matrix",a,ndim,n);
/* showv("solve: RHS",b,n); /**/
/* forward elimination */
if(n!=1)
{nm1=n-1;
for (k=0; k {kp1=k+1;
m=ipvt[k];
t=b[m];
b[m]=b[k];
b[k]=t;
for (i=kp1; i b[i] += aa(i,k)*t;
}
/* printf("b\n");
/* printf("n=%d, ipvt= ", n); for (k=0; k /* printf("\n");
/* for (k=0; k /* {for (i=0; i /* printf("%f ", aa(k, i));
/* printf("\n");
/* }
/*
/* showv("\nafter forward elimination",b,n); /**/
/* back substitution */
for (k=nm1; k; k--)
{b[k] /= aa(k,k);
t=-b[k];
for (i=0; i b[i] += aa(i,k)*t;
}
}
b[0] /= a[0];
/* printf("c");
/* showv("\nafter back substitution",b,n); /**/
/*
free(ix);
*/
}

#ifdef TEST

/* sample program for decomp and solve */

main()
{ double x[13][13], p[13][13], a[13][13], b[13], cond, condp1;
int ipvt[13], i, j, k, n, ndim;

ndim=13;
n=3;
a[0][0]=10.;
a[1][0]=-3.;
a[2][0]= 5.;
a[0][1]=-7.;
a[1][1]= 2.;
a[2][1]=-1.;
a[0][2]= 0.;
a[1][2]= 6.;
a[2][2]= 5.;
for (i=0; i {for (j=0; j printf("%8.4f ",a[i][j]);
printf("\n");
}
printf("\n");
cond=decomp(ndim,n,a,ipvt);
printf("cond = %f\n\n",cond);
printf("\nipvt\n");
for (i=0; i printf("%8d \n",ipvt[i]);
condp1=cond+1.;
if(condp1==cond) exit();
b[0]=7.;
b[1]=4.;
b[2]=6.;
showv("RHS",b,n);
solve(ndim,n,a,b,ipvt);
showv("solution",b,n);
a[0][0]=10.;
a[1][0]=-3.;
a[2][0]= 5.;
a[0][1]=-7.;
a[1][1]= 2.;
a[2][1]=-1.;
a[0][2]= 0.;
a[1][2]= 6.;
a[2][2]= 5.;
showm("a, matrix to be inverted",a,ndim,n);
cond=invert(ndim,n,a,x);
printf("cond = %f\n\n",cond);
printf("ipvt\n");
for (i=0; i printf("%8d \n",ipvt[i]);
condp1=cond+1.;
if(condp1==cond) exit();
showm("x, inverse",x,ndim,n);
a[0][0]=10.;
a[1][0]=-3.;
a[2][0]= 5.;
a[0][1]=-7.;
a[1][1]= 2.;
a[2][1]=-1.;
a[0][2]= 0.;
a[1][2]= 6.;
a[2][2]= 5.;
for (k=0; k for (j=0; j {p[k][j]=0.;
for (i=0; i p[k][j] += a[k][i]*x[i][j];
}
showm("a*x",p,ndim,n);
for (k=0; k for (j=0; j {p[k][j]=0.;
for (i=0; i p[k][j] += x[k][i]*a[i][j];
}
showm("x*a",p,ndim,n);
}

showv(s,a,n) char *s; double a[]; int n;
{ int i,j;
printf("%s at %04xH\n",s,a);
for (i=0; i printf("%8.4f \n",a[i]);
printf("\n");
}

showiv(s,a,n) char *s; int a[]; int n;
{ int i,j,fail=0;
printf("%s at %04xH",s,a);
for (i=0; i {printf("%8d ",a[i]);
if(a[i]>=n) fail=1;
}
printf("\n");
if(fail) exit();
}

showm(s,a,ndim,n) char *s; double a[]; int ndim,n;
{ int i,j;
printf("%s\n",s);
for (i=0; i {for (j=0; j printf("%8.4f ",a[i*ndim+j]);
printf("\n");
}
printf("\n");
}
#endif


  3 Responses to “Category : Science and Education
Archive   : LAPLACE.ZIP
Filename : SOLVE.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/