1 #ifndef lint 2 static char vcid[] = "$Id: dgefa3.c,v 1.2 1996/04/28 02:24:51 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Inverts 3 by 3 matrix using partial pivoting. 6 */ 7 #include "petsc.h" 8 9 int Kernel_A_gets_inverse_A_3(Scalar *a) 10 { 11 int i__2, i__3, kp1, j, k, l,ll,i,ipvt_l[3],*ipvt = ipvt_l,kb,k3; 12 int k4,j3; 13 Scalar t,*aa,*ax,*ay,work_l[9],*work = work_l,stmp; 14 double tmp,max; 15 16 /* gaussian elimination with partial pivoting */ 17 18 /* Parameter adjustments */ 19 --ipvt; 20 a -= 4; 21 22 for (k = 1; k <= 2; ++k) { 23 kp1 = k + 1; 24 k3 = 3*k; 25 k4 = k3 + k; 26 /* find l = pivot index */ 27 28 i__2 = 4 - k; 29 aa = &a[k4]; 30 max = PetscAbsScalar(aa[0]); 31 l = 1; 32 for ( ll=1; ll<i__2; ll++ ) { 33 tmp = PetscAbsScalar(aa[ll]); 34 if (tmp > max) { max = tmp; l = ll+1;} 35 } 36 l += k - 1; 37 ipvt[k] = l; 38 39 if (a[l + k3] == 0.) { 40 SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 41 } 42 43 /* interchange if necessary */ 44 45 if (l != k) { 46 t = a[l + k3]; 47 a[l + k3] = a[k4]; 48 a[k4] = t; 49 } 50 51 /* compute multipliers */ 52 53 t = -1. / a[k4]; 54 i__2 = 3 - k; 55 aa = &a[1 + k4]; 56 for ( ll=0; ll<i__2; ll++ ) { 57 aa[ll] *= t; 58 } 59 60 /* row elimination with column indexing */ 61 62 ax = &a[k4+1]; 63 for (j = kp1; j <= 3; ++j) { 64 j3 = 3*j; 65 t = a[l + j3]; 66 if (l != k) { 67 a[l + j3] = a[k + j3]; 68 a[k + j3] = t; 69 } 70 71 i__3 = 3 - k; 72 ay = &a[k+1+j3]; 73 for ( ll=0; ll<i__3; ll++ ) { 74 ay[ll] += t*ax[ll]; 75 } 76 } 77 } 78 ipvt[3] = 3; 79 if (a[12] == 0.) { 80 SETERRQ(3,"Linpack_DGEFA:Zero pivot,final row"); 81 } 82 83 /* 84 Now form the inverse 85 */ 86 87 --work; 88 89 /* compute inverse(u) */ 90 91 for (k = 1; k <= 3; ++k) { 92 k3 = 3*k; 93 k4 = k3 + k; 94 a[k4] = 1.0 / a[k4]; 95 t = -a[k4]; 96 i__2 = k - 1; 97 aa = &a[k3 + 1]; 98 for ( ll=0; ll<i__2; ll++ ) aa[ll] *= t; 99 kp1 = k + 1; 100 if (3 < kp1) continue; 101 ax = aa; 102 for (j = kp1; j <= 3; ++j) { 103 j3 = 3*j; 104 t = a[k + j3]; 105 a[k + j3] = 0.; 106 ay = &a[j3 + 1]; 107 for ( ll=0; ll<k; ll++ ) { 108 ay[ll] += t*ax[ll]; 109 } 110 } 111 } 112 113 /* form inverse(u)*inverse(l) */ 114 115 for (kb = 1; kb <= 2; ++kb) { 116 k = 3 - kb; 117 k3 = 3*k; 118 kp1 = k + 1; 119 aa = a + k3; 120 for (i = kp1; i <= 3; ++i) { 121 work[i] = aa[i]; 122 aa[i] = 0.0; 123 } 124 for (j = kp1; j <= 3; ++j) { 125 t = work[j]; 126 ax = &a[3*j + 1]; 127 ay = &a[k3 + 1]; 128 ay[0] += t*ax[0]; 129 ay[1] += t*ax[1]; 130 ay[2] += t*ax[2]; 131 } 132 l = ipvt[k]; 133 if (l != k) { 134 ax = &a[k3 + 1]; 135 ay = &a[3*l + 1]; 136 stmp = ax[0]; ax[0] = ay[0]; ay[0] = stmp; 137 stmp = ax[1]; ax[1] = ay[1]; ay[1] = stmp; 138 stmp = ax[2]; ax[2] = ay[2]; ay[2] = stmp; 139 } 140 } 141 return 0; 142 } 143 144