1*2bb67c89SBarry Smith /* dgedi.f -- translated by f2c (version of 25 March 1992 12:58:56). 2*2bb67c89SBarry Smith You must link the resulting object file with the libraries: 3*2bb67c89SBarry Smith -lF77 -lI77 -lm -lc (in that order) 4*2bb67c89SBarry Smith */ 5*2bb67c89SBarry Smith 6*2bb67c89SBarry Smith #include <f2c.h> 7*2bb67c89SBarry Smith 8*2bb67c89SBarry Smith /* Table of constant values */ 9*2bb67c89SBarry Smith 10*2bb67c89SBarry Smith static integer c__1 = 1; 11*2bb67c89SBarry Smith 12*2bb67c89SBarry Smith /* Subroutine */ int dgedi_(a, lda, n, ipvt, det, work, job) 13*2bb67c89SBarry Smith doublereal *a; 14*2bb67c89SBarry Smith integer *lda, *n, *ipvt; 15*2bb67c89SBarry Smith doublereal *det, *work; 16*2bb67c89SBarry Smith integer *job; 17*2bb67c89SBarry Smith { 18*2bb67c89SBarry Smith /* System generated locals */ 19*2bb67c89SBarry Smith integer a_dim1, a_offset, i__1, i__2; 20*2bb67c89SBarry Smith 21*2bb67c89SBarry Smith /* Local variables */ 22*2bb67c89SBarry Smith static integer i, j, k, l; 23*2bb67c89SBarry Smith static doublereal t; 24*2bb67c89SBarry Smith extern /* Subroutine */ int dscal_(), dswap_(), daxpy_(); 25*2bb67c89SBarry Smith static integer kb, kp1, nm1; 26*2bb67c89SBarry Smith 27*2bb67c89SBarry Smith 28*2bb67c89SBarry Smith /* dgedi computes the determinant and inverse of a matrix */ 29*2bb67c89SBarry Smith /* using the factors computed by dgeco or dgefa. */ 30*2bb67c89SBarry Smith 31*2bb67c89SBarry Smith /* on entry */ 32*2bb67c89SBarry Smith 33*2bb67c89SBarry Smith /* a double precision(lda, n) */ 34*2bb67c89SBarry Smith /* the output from dgeco or dgefa. */ 35*2bb67c89SBarry Smith 36*2bb67c89SBarry Smith /* lda integer */ 37*2bb67c89SBarry Smith /* the leading dimension of the array a . */ 38*2bb67c89SBarry Smith 39*2bb67c89SBarry Smith /* n integer */ 40*2bb67c89SBarry Smith /* the order of the matrix a . */ 41*2bb67c89SBarry Smith 42*2bb67c89SBarry Smith /* ipvt integer(n) */ 43*2bb67c89SBarry Smith /* the pivot vector from dgeco or dgefa. */ 44*2bb67c89SBarry Smith 45*2bb67c89SBarry Smith /* work double precision(n) */ 46*2bb67c89SBarry Smith /* work vector. contents destroyed. */ 47*2bb67c89SBarry Smith 48*2bb67c89SBarry Smith /* job integer */ 49*2bb67c89SBarry Smith /* = 11 both determinant and inverse. */ 50*2bb67c89SBarry Smith /* = 01 inverse only. */ 51*2bb67c89SBarry Smith /* = 10 determinant only. */ 52*2bb67c89SBarry Smith 53*2bb67c89SBarry Smith /* on return */ 54*2bb67c89SBarry Smith 55*2bb67c89SBarry Smith /* a inverse of original matrix if requested. */ 56*2bb67c89SBarry Smith /* otherwise unchanged. */ 57*2bb67c89SBarry Smith 58*2bb67c89SBarry Smith /* det double precision(2) */ 59*2bb67c89SBarry Smith /* determinant of original matrix if requested. */ 60*2bb67c89SBarry Smith /* otherwise not referenced. */ 61*2bb67c89SBarry Smith /* determinant = det(1) * 10.0**det(2) */ 62*2bb67c89SBarry Smith /* with 1.0 .le. dabs(det(1)) .lt. 10.0 */ 63*2bb67c89SBarry Smith /* or det(1) .eq. 0.0 . */ 64*2bb67c89SBarry Smith 65*2bb67c89SBarry Smith /* error condition */ 66*2bb67c89SBarry Smith 67*2bb67c89SBarry Smith /* a division by zero will occur if the input factor contains */ 68*2bb67c89SBarry Smith /* a zero on the diagonal and the inverse is requested. */ 69*2bb67c89SBarry Smith /* it will not occur if the subroutines are called correctly */ 70*2bb67c89SBarry Smith /* and if dgeco has set rcond .gt. 0.0 or dgefa has set */ 71*2bb67c89SBarry Smith /* info .eq. 0 . */ 72*2bb67c89SBarry Smith 73*2bb67c89SBarry Smith /* linpack. this version dated 08/14/78 . */ 74*2bb67c89SBarry Smith /* cleve moler, university of new mexico, argonne national lab. */ 75*2bb67c89SBarry Smith 76*2bb67c89SBarry Smith /* subroutines and functions */ 77*2bb67c89SBarry Smith 78*2bb67c89SBarry Smith /* blas daxpy,dscal,dswap */ 79*2bb67c89SBarry Smith /* fortran dabs,mod */ 80*2bb67c89SBarry Smith 81*2bb67c89SBarry Smith /* internal variables */ 82*2bb67c89SBarry Smith 83*2bb67c89SBarry Smith 84*2bb67c89SBarry Smith 85*2bb67c89SBarry Smith /* compute inverse(u) */ 86*2bb67c89SBarry Smith 87*2bb67c89SBarry Smith /* Parameter adjustments */ 88*2bb67c89SBarry Smith --work; 89*2bb67c89SBarry Smith --det; 90*2bb67c89SBarry Smith --ipvt; 91*2bb67c89SBarry Smith a_dim1 = *lda; 92*2bb67c89SBarry Smith a_offset = a_dim1 + 1; 93*2bb67c89SBarry Smith a -= a_offset; 94*2bb67c89SBarry Smith 95*2bb67c89SBarry Smith /* Function Body */ 96*2bb67c89SBarry Smith if (*job % 10 == 0) { 97*2bb67c89SBarry Smith goto L150; 98*2bb67c89SBarry Smith } 99*2bb67c89SBarry Smith i__1 = *n; 100*2bb67c89SBarry Smith for (k = 1; k <= i__1; ++k) { 101*2bb67c89SBarry Smith a[k + k * a_dim1] = 1. / a[k + k * a_dim1]; 102*2bb67c89SBarry Smith t = -a[k + k * a_dim1]; 103*2bb67c89SBarry Smith i__2 = k - 1; 104*2bb67c89SBarry Smith dscal_(&i__2, &t, &a[k * a_dim1 + 1], &c__1); 105*2bb67c89SBarry Smith kp1 = k + 1; 106*2bb67c89SBarry Smith if (*n < kp1) { 107*2bb67c89SBarry Smith goto L90; 108*2bb67c89SBarry Smith } 109*2bb67c89SBarry Smith i__2 = *n; 110*2bb67c89SBarry Smith for (j = kp1; j <= i__2; ++j) { 111*2bb67c89SBarry Smith t = a[k + j * a_dim1]; 112*2bb67c89SBarry Smith a[k + j * a_dim1] = 0.; 113*2bb67c89SBarry Smith daxpy_(&k, &t, &a[k * a_dim1 + 1], &c__1, &a[j * a_dim1 + 1], & 114*2bb67c89SBarry Smith c__1); 115*2bb67c89SBarry Smith /* L80: */ 116*2bb67c89SBarry Smith } 117*2bb67c89SBarry Smith L90: 118*2bb67c89SBarry Smith /* L100: */ 119*2bb67c89SBarry Smith ; 120*2bb67c89SBarry Smith } 121*2bb67c89SBarry Smith 122*2bb67c89SBarry Smith /* form inverse(u)*inverse(l) */ 123*2bb67c89SBarry Smith 124*2bb67c89SBarry Smith nm1 = *n - 1; 125*2bb67c89SBarry Smith if (nm1 < 1) { 126*2bb67c89SBarry Smith goto L140; 127*2bb67c89SBarry Smith } 128*2bb67c89SBarry Smith i__1 = nm1; 129*2bb67c89SBarry Smith for (kb = 1; kb <= i__1; ++kb) { 130*2bb67c89SBarry Smith k = *n - kb; 131*2bb67c89SBarry Smith kp1 = k + 1; 132*2bb67c89SBarry Smith i__2 = *n; 133*2bb67c89SBarry Smith for (i = kp1; i <= i__2; ++i) { 134*2bb67c89SBarry Smith work[i] = a[i + k * a_dim1]; 135*2bb67c89SBarry Smith a[i + k * a_dim1] = 0.; 136*2bb67c89SBarry Smith /* L110: */ 137*2bb67c89SBarry Smith } 138*2bb67c89SBarry Smith i__2 = *n; 139*2bb67c89SBarry Smith for (j = kp1; j <= i__2; ++j) { 140*2bb67c89SBarry Smith t = work[j]; 141*2bb67c89SBarry Smith daxpy_(n, &t, &a[j * a_dim1 + 1], &c__1, &a[k * a_dim1 + 1], & 142*2bb67c89SBarry Smith c__1); 143*2bb67c89SBarry Smith /* L120: */ 144*2bb67c89SBarry Smith } 145*2bb67c89SBarry Smith l = ipvt[k]; 146*2bb67c89SBarry Smith if (l != k) { 147*2bb67c89SBarry Smith dswap_(n, &a[k * a_dim1 + 1], &c__1, &a[l * a_dim1 + 1], &c__1); 148*2bb67c89SBarry Smith } 149*2bb67c89SBarry Smith /* L130: */ 150*2bb67c89SBarry Smith } 151*2bb67c89SBarry Smith L140: 152*2bb67c89SBarry Smith L150: 153*2bb67c89SBarry Smith return 0; 154*2bb67c89SBarry Smith } /* dgedi_ */ 155*2bb67c89SBarry Smith 156