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