1 2 /* 3 This routine was converted by f2c from Linpack source 4 linpack. this version dated 08/14/78 5 cleve moler, university of new mexico, argonne national lab. 6 */ 7 #include "petsc.h" 8 9 int Linpack_DGEFA(Scalar *a, int n, int *ipvt) 10 { 11 int i__2, i__3, kp1, nm1, j, k, l,ll; 12 Scalar t,*aa,*ax,*ay; 13 double tmp,max; 14 15 /* gaussian elimination with partial pivoting */ 16 17 /* Parameter adjustments */ 18 --ipvt; 19 a -= n + 1; 20 21 /* Function Body */ 22 nm1 = n - 1; 23 if (nm1 < 1) { 24 goto L70; 25 } 26 for (k = 1; k <= nm1; ++k) { 27 kp1 = k + 1; 28 29 /* find l = pivot index */ 30 31 i__2 = n - k + 1; 32 aa = &a[k + k * n]; 33 max = PetscAbsScalar(aa[0]); 34 l = 1; 35 for ( ll=1; ll<i__2; ll++ ) { 36 tmp = PetscAbsScalar(aa[ll]); 37 if (tmp > max) { max = tmp; l = ll+1;} 38 } 39 l += k - 1; 40 ipvt[k] = l; 41 42 if (a[l + k * n] == 0.) { 43 SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 44 } 45 46 /* interchange if necessary */ 47 48 if (l != k) { 49 t = a[l + k * n]; 50 a[l + k * n] = a[k + k * n]; 51 a[k + k * n] = t; 52 } 53 54 /* compute multipliers */ 55 56 t = -1. / a[k + k * n]; 57 i__2 = n - k; 58 aa = &a[k + 1 + k * n]; 59 for ( ll=0; ll<i__2; ll++ ) { 60 aa[ll] *= t; 61 } 62 63 /* row elimination with column indexing */ 64 65 ax = &a[k+1+k*n]; 66 for (j = kp1; j <= n; ++j) { 67 t = a[l + j * n]; 68 if (l != k) { 69 a[l + j * n] = a[k + j * n]; 70 a[k + j * n] = t; 71 } 72 73 i__3 = n - k; 74 ay = &a[k+1+j*n]; 75 for ( ll=0; ll<i__3; ll++ ) { 76 ay[ll] += t*ax[ll]; 77 } 78 } 79 } 80 L70: 81 ipvt[n] = n; 82 if (a[n + n * n] == 0.) { 83 SETERRQ(n,"Linpack_DGEFA:Zero pivot,final row"); 84 } 85 return 0; 86 } 87 88 /* 89 Version for 3 by 3 matrices. 90 91 THIS CODE HAS NOT BEEN TESTED YET! 92 93 */ 94 int Linpack_DGEFA_3(Scalar *a, int *ipvt) 95 { 96 int i__2, i__3, kp1, j, k, l,ll; 97 Scalar t,*aa,*ax,*ay; 98 double tmp,max; 99 100 /* gaussian elimination with partial pivoting */ 101 102 /* Parameter adjustments */ 103 --ipvt; 104 a -= 4; 105 106 /* Function Body */ 107 for (k = 1; k <= 2; ++k) { 108 kp1 = k + 1; 109 110 /* find l = pivot index */ 111 112 i__2 = 4 - k; 113 aa = &a[4*k]; 114 max = PetscAbsScalar(aa[0]); 115 l = 1; 116 for ( ll=1; ll<i__2; ll++ ) { 117 tmp = PetscAbsScalar(aa[ll]); 118 if (tmp > max) { max = tmp; l = ll+1;} 119 } 120 l += k - 1; 121 ipvt[k] = l; 122 123 if (a[l + 3*k] == 0.) { 124 SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 125 } 126 127 /* interchange if necessary */ 128 129 if (l != k) { 130 t = a[l + 3*k]; 131 a[l + 3*k] = a[4*k]; 132 a[4*k] = t; 133 } 134 135 /* compute multipliers */ 136 137 t = -1. / a[4*k]; 138 i__2 = 3 - k; 139 aa = &a[1 + 4*k]; 140 for ( ll=0; ll<i__2; ll++ ) { 141 aa[ll] *= t; 142 } 143 144 /* row elimination with column indexing */ 145 146 ax = &a[4*k+1]; 147 for (j = kp1; j <= 3; ++j) { 148 t = a[l + 3*j]; 149 if (l != k) { 150 a[l + 3*j] = a[k + 3*j]; 151 a[k + 3*j] = t; 152 } 153 154 i__3 = 3 - k; 155 ay = &a[k+1+3*j]; 156 for ( ll=0; ll<i__3; ll++ ) { 157 ay[ll] += t*ax[ll]; 158 } 159 } 160 } 161 ipvt[3] = 3; 162 if (a[12] == 0.) { 163 SETERRQ(3,"Linpack_DGEFA:Zero pivot,final row"); 164 } 165 return 0; 166 } 167 168 169 #ifndef lint 170 static char vcid[] = "$Id: dgedi.c,v 1.2 1996/03/04 05:16:18 bsmith Exp bsmith $"; 171 #endif 172 173 /* 174 This file creating by running f2c 175 linpack. this version dated 08/14/78 176 cleve moler, university of new mexico, argonne national lab. 177 178 Computes the inverse of a matrix given its factors and pivots 179 calculated by Linpack_DGEFA(). 180 */ 181 182 #include "petsc.h" 183 184 int Linpack_DGEDI(Scalar *a,int n,int *ipvt,Scalar *work) 185 { 186 int i__2,kb, kp1, nm1,i, j, k, l, ll; 187 Scalar t, *aa,*ax,*ay,tmp; 188 189 --work; 190 --ipvt; 191 a -= n + 1; 192 193 /* compute inverse(u) */ 194 195 for (k = 1; k <= n; ++k) { 196 a[k + k * n] = 1. / a[k + k * n]; 197 t = -a[k + k * n]; 198 i__2 = k - 1; 199 aa = &a[k * n + 1]; 200 for ( ll=0; ll<i__2; ll++ ) aa[ll] *= t; 201 kp1 = k + 1; 202 if (n < kp1) continue; 203 ax = aa; 204 for (j = kp1; j <= n; ++j) { 205 t = a[k + j * n]; 206 a[k + j * n] = 0.; 207 ay = &a[j * n + 1]; 208 for ( ll=0; ll<k; ll++ ) { 209 ay[ll] += t*ax[ll]; 210 } 211 } 212 } 213 214 /* form inverse(u)*inverse(l) */ 215 216 nm1 = n - 1; 217 if (nm1 < 1) { 218 return 0; 219 } 220 for (kb = 1; kb <= nm1; ++kb) { 221 k = n - kb; 222 kp1 = k + 1; 223 aa = a + k * n; 224 for (i = kp1; i <= n; ++i) { 225 work[i] = aa[i]; 226 aa[i] = 0.; 227 } 228 for (j = kp1; j <= n; ++j) { 229 t = work[j]; 230 ax = &a[j * n + 1]; 231 ay = &a[k * n + 1]; 232 for ( ll=0; ll<n; ll++ ) { 233 ay[ll] += t*ax[ll]; 234 } 235 } 236 l = ipvt[k]; 237 if (l != k) { 238 ax = &a[k * n + 1]; 239 ay = &a[l * n + 1]; 240 for ( ll=0; ll<n; ll++ ) { 241 tmp = ax[ll]; 242 ax[ll] = ay[ll]; 243 ay[ll] = tmp; 244 } 245 } 246 } 247 return 0; 248 } 249 250