| 1 | #include "slalib.h"
|
|---|
| 2 | #include "slamac.h"
|
|---|
| 3 | void slaDmat ( int n, double *a, double *y, double *d, int *jf, int *iw)
|
|---|
| 4 | /*
|
|---|
| 5 | ** - - - - - - - -
|
|---|
| 6 | ** s l a D m a t
|
|---|
| 7 | ** - - - - - - - -
|
|---|
| 8 | **
|
|---|
| 9 | ** Matrix inversion & solution of simultaneous equations.
|
|---|
| 10 | **
|
|---|
| 11 | ** (double precision)
|
|---|
| 12 | **
|
|---|
| 13 | ** For the set of n simultaneous equations in n unknowns:
|
|---|
| 14 | ** a.y = x
|
|---|
| 15 | **
|
|---|
| 16 | ** where:
|
|---|
| 17 | ** a is a non-singular n x n matrix
|
|---|
| 18 | ** y is the vector of n unknowns
|
|---|
| 19 | ** x is the known vector
|
|---|
| 20 | **
|
|---|
| 21 | ** slaDmat computes:
|
|---|
| 22 | ** the inverse of matrix a
|
|---|
| 23 | ** the determinant of matrix a
|
|---|
| 24 | ** the vector of n unknowns
|
|---|
| 25 | **
|
|---|
| 26 | ** Arguments:
|
|---|
| 27 | **
|
|---|
| 28 | ** symbol type dimension before after
|
|---|
| 29 | **
|
|---|
| 30 | ** n int no. of unknowns unchanged
|
|---|
| 31 | ** *a double [n][n] matrix inverse
|
|---|
| 32 | ** *y double [n] vector solution
|
|---|
| 33 | ** *d double - determinant
|
|---|
| 34 | ** > *jf int - singularity flag
|
|---|
| 35 | ** *iw int [n] - workspace
|
|---|
| 36 | **
|
|---|
| 37 | **
|
|---|
| 38 | ** > jf is the singularity flag. If the matrix is non-singular,
|
|---|
| 39 | ** jf=0 is returned. If the matrix is singular, jf=-1 & d=0.0 are
|
|---|
| 40 | ** returned. In the latter case, the contents of array a on return
|
|---|
| 41 | ** are undefined.
|
|---|
| 42 | **
|
|---|
| 43 | ** Algorithm:
|
|---|
| 44 | ** Gaussian elimination with partial pivoting.
|
|---|
| 45 | **
|
|---|
| 46 | ** Speed:
|
|---|
| 47 | ** Very fast.
|
|---|
| 48 | **
|
|---|
| 49 | ** Accuracy:
|
|---|
| 50 | ** Fairly accurate - errors 1 to 4 times those of routines optimized
|
|---|
| 51 | ** for accuracy.
|
|---|
| 52 | **
|
|---|
| 53 | ** Example call (note handling of "adjustable dimension" 2D array):
|
|---|
| 54 | **
|
|---|
| 55 | ** double a[MP][MP], v[MP], d;
|
|---|
| 56 | ** int j, iw[MP];
|
|---|
| 57 | ** :
|
|---|
| 58 | ** slaDmat ( n, (double *) a, v, &d, &j, iw );
|
|---|
| 59 | **
|
|---|
| 60 | ** Last revision: 20 February 1995
|
|---|
| 61 | **
|
|---|
| 62 | ** Copyright P.T.Wallace. All rights reserved.
|
|---|
| 63 | */
|
|---|
| 64 |
|
|---|
| 65 | #define TINY 1e-20
|
|---|
| 66 |
|
|---|
| 67 | {
|
|---|
| 68 | int k, imx, i, j, ki;
|
|---|
| 69 | double amx, t, yk;
|
|---|
| 70 |
|
|---|
| 71 | /* Pointers to beginnings of rows in matrix a[n][n] */
|
|---|
| 72 |
|
|---|
| 73 | double *ak, /* row k */
|
|---|
| 74 | *ai, /* row i */
|
|---|
| 75 | *aimx; /* row imx */
|
|---|
| 76 |
|
|---|
| 77 | *jf = 0;
|
|---|
| 78 | *d = 1.0;
|
|---|
| 79 | for ( k = 0, ak = a; k < n; k++, ak += n ) {
|
|---|
| 80 | amx = fabs ( ak[k] );
|
|---|
| 81 | imx = k;
|
|---|
| 82 | aimx = ak;
|
|---|
| 83 | if ( k != n ) {
|
|---|
| 84 | for ( i = k + 1, ai = ak + n; i < n; i++, ai += n ) {
|
|---|
| 85 | t = fabs ( ai[k] );
|
|---|
| 86 | if ( t > amx ) {
|
|---|
| 87 | amx = t;
|
|---|
| 88 | imx = i;
|
|---|
| 89 | aimx = ai;
|
|---|
| 90 | }
|
|---|
| 91 | }
|
|---|
| 92 | }
|
|---|
| 93 | if ( amx < TINY ) {
|
|---|
| 94 | *jf = -1;
|
|---|
| 95 | } else {
|
|---|
| 96 | if ( imx != k ) {
|
|---|
| 97 | for ( j = 0; j < n; j++ ) {
|
|---|
| 98 | t = ak[j];
|
|---|
| 99 | ak[j] = aimx[j];
|
|---|
| 100 | aimx[j] = t;
|
|---|
| 101 | }
|
|---|
| 102 | t = y[k];
|
|---|
| 103 | y[k] = y[imx];
|
|---|
| 104 | y[imx] = t;
|
|---|
| 105 | *d = - *d;
|
|---|
| 106 | }
|
|---|
| 107 | iw[k] = imx;
|
|---|
| 108 | *d *= ak[k];
|
|---|
| 109 | if ( fabs ( *d ) < TINY ) {
|
|---|
| 110 | *jf = -1;
|
|---|
| 111 | } else {
|
|---|
| 112 | ak[k] = 1.0 / ak[k];
|
|---|
| 113 | for ( j = 0; j < n; j++ ) {
|
|---|
| 114 | if ( j != k ) {
|
|---|
| 115 | ak[j] *= ak[k];
|
|---|
| 116 | }
|
|---|
| 117 | }
|
|---|
| 118 | yk = y[k] * ak[k];
|
|---|
| 119 | y[k] = yk;
|
|---|
| 120 | for ( i = 0, ai = a; i < n; i++, ai += n ) {
|
|---|
| 121 | if ( i != k ) {
|
|---|
| 122 | for ( j = 0; j < n; j++ ) {
|
|---|
| 123 | if ( j != k ) {
|
|---|
| 124 | ai[j] -= ai[k] * ak[j];
|
|---|
| 125 | }
|
|---|
| 126 | }
|
|---|
| 127 | y[i] -= ai[k] * yk;
|
|---|
| 128 | }
|
|---|
| 129 | }
|
|---|
| 130 | for ( i = 0, ai = a; i < n; i++, ai += n ) {
|
|---|
| 131 | if ( i != k ) {
|
|---|
| 132 | ai[k] *= - ak[k];
|
|---|
| 133 | }
|
|---|
| 134 | }
|
|---|
| 135 | }
|
|---|
| 136 | }
|
|---|
| 137 | }
|
|---|
| 138 | if ( *jf != 0 ) {
|
|---|
| 139 | *d = 0.0;
|
|---|
| 140 | } else {
|
|---|
| 141 | for ( k = n; k-- > 0; ) {
|
|---|
| 142 | ki = iw[k];
|
|---|
| 143 | if ( k != ki ) {
|
|---|
| 144 | for ( i = 0, ai = a; i < n; i++, ai += n ) {
|
|---|
| 145 | t = ai[k];
|
|---|
| 146 | ai[k] = ai[ki];
|
|---|
| 147 | ai[ki] = t;
|
|---|
| 148 | }
|
|---|
| 149 | }
|
|---|
| 150 | }
|
|---|
| 151 | }
|
|---|
| 152 | }
|
|---|
