/*****************************************************************************/ /* */ /* spline_miv.c */ /* */ /* "spline_miv" inverts a square matrix using Gauss-Jordan elimination. */ /* This routine is based on a VAX/VMS assembly code version which in turn */ /* had many previous incarnations. */ /* */ /* Arguments: n dimension of square matrices (a[n,n], b[n,n]) */ /* a matrix to be inverted. 'a' is destroyed by routine. */ /* b inverse matrix */ /* m temporary storage vector, dimension is m[n] */ /* */ /* Return: error code. If negative, a computation error was */ /* encountered. */ /* */ /* Carl W. Akerlof */ /* Center for Particle Astrophysics */ /* 301 Le Conte Hall */ /* University of California */ /* Berkeley, California 94720 */ /* */ /* June 29, 1993 */ /* */ /*****************************************************************************/ #define ABS(A) ((A >= 0.0) ? (A) : (-(A))) int spline_miv(int n, float a[], float b[], int m[]) { int i, ij, j, k, nn, p; float u, v; nn=n*n; /*****************************************************************************/ /* */ /* Initialize matrix B */ /* */ /*****************************************************************************/ for (i=0; i < nn; i++) { b[i]=0.0; } for (i=0; i < n; i++) { b[(n+1)*i]=1.0; m[i]=i; } /*****************************************************************************/ /* */ /* Main matrix inversion processing loop */ /* */ /*****************************************************************************/ for (p=n; p > 0; p--) { /*****************************************************************************/ /* */ /* Look for largest element in unreduced columns of matrix A */ /* */ /*****************************************************************************/ u=0.0; for (i=0; i < p; i++) { for (j=n*m[i]; j < n*(m[i]+1); j++) { if (ABS(a[j]) > u) { k=j; u=ABS(a[j]); } } } if (u <= 0.0) { return (-1); } i=k%n; j=k-i; i*=n; u=a[k]; /*****************************************************************************/ /* */ /* Exchange columns I and J of A */ /* Divide column I of A by A(I,J) */ /* Exchange columns I and J of B */ /* Divide column I of B by A(I,J) */ /* */ /*****************************************************************************/ for (k=0; k < n; k++) { v=a[j]; a[j]=a[i]; a[i]=v/u; v=b[j]; b[j++]=b[i]; b[i++]=v/u; } i-=n; j-=n; /*****************************************************************************/ /* */ /* Find and replace element of M( ) to be reduced */ /* */ /*****************************************************************************/ for (k=0; k < p; k++) { if (m[k] == i/n) { break; } } m[k]=m[p-1]; /*****************************************************************************/ /* */ /* Clear I'th row of A( , ) by subtracting column I from column J for */ /* all J (I != J) */ /* Perform identical operations on matrix B */ /* */ /*****************************************************************************/ for (j=0; j < nn;) { if (j != i) { ij=i/n+j; u=a[ij]; a[ij]=0.0; for (k=0; k < n; k++) { if (j != ij) { a[j]-=u*a[i]; } b[j++]-=u*b[i++]; } i-=n; } else { j+=n; } } } return (0); }