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