1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dgefa2.c,v 1.1 1999/03/16 22:23:56 bsmith Exp balay $"; 3 #endif 4 /* 5 Inverts 2 by 2 matrix using partial pivoting. 6 7 Used by the sparse factorization routines in 8 src/mat/impls/baij/seq and src/mat/impls/bdiag/seq 9 10 See also src/inline/ilu.h 11 12 This is a combination of the Linpack routines 13 dgefa() and dgedi() specialized for a size of 2. 14 15 */ 16 #include "petsc.h" 17 18 #undef __FUNC__ 19 #define __FUNC__ "Kernel_A_gets_inverse_A_2" 20 int Kernel_A_gets_inverse_A_2(MatScalar *a) 21 { 22 int i__2, i__3, kp1, j, k, l,ll,i,ipvt_l[2],*ipvt = ipvt_l-1,k3; 23 int k4,j3; 24 MatScalar *aa,*ax,*ay,work_l[4],*work = work_l-1,stmp; 25 MatFloat tmp,max; 26 27 /* gaussian elimination with partial pivoting */ 28 29 PetscFunctionBegin; 30 /* Parameter adjustments */ 31 a -= 3; 32 33 /*for (k = 1; k <= 1; ++k) {*/ 34 k = 1; 35 kp1 = k + 1; 36 k3 = 2*k; 37 k4 = k3 + k; 38 /* find l = pivot index */ 39 40 i__2 = 2 - k; 41 aa = &a[k4]; 42 max = PetscAbsScalar(aa[0]); 43 l = 1; 44 for ( ll=1; ll<i__2; ll++ ) { 45 tmp = PetscAbsScalar(aa[ll]); 46 if (tmp > max) { max = tmp; l = ll+1;} 47 } 48 l += k - 1; 49 ipvt[k] = l; 50 51 if (a[l + k3] == 0.) { 52 SETERRQ(k,0,"Zero pivot"); 53 } 54 55 /* interchange if necessary */ 56 57 if (l != k) { 58 stmp = a[l + k3]; 59 a[l + k3] = a[k4]; 60 a[k4] = stmp; 61 } 62 63 /* compute multipliers */ 64 65 stmp = -1. / a[k4]; 66 i__2 = 2 - k; 67 aa = &a[1 + k4]; 68 for ( ll=0; ll<i__2; ll++ ) { 69 aa[ll] *= stmp; 70 } 71 72 /* row elimination with column indexing */ 73 74 ax = &a[k4+1]; 75 for (j = kp1; j <= 2; ++j) { 76 j3 = 2*j; 77 stmp = a[l + j3]; 78 if (l != k) { 79 a[l + j3] = a[k + j3]; 80 a[k + j3] = stmp; 81 } 82 83 i__3 = 2 - k; 84 ay = &a[1+k+j3]; 85 for ( ll=0; ll<i__3; ll++ ) { 86 ay[ll] += stmp*ax[ll]; 87 } 88 } 89 /*}*/ 90 ipvt[2] = 2; 91 if (a[6] == 0.) { 92 SETERRQ(3,0,"Zero pivot,final row"); 93 } 94 95 /* 96 Now form the inverse 97 */ 98 99 /* compute inverse(u) */ 100 101 for (k = 1; k <= 2; ++k) { 102 k3 = 2*k; 103 k4 = k3 + k; 104 a[k4] = 1.0 / a[k4]; 105 stmp = -a[k4]; 106 i__2 = k - 1; 107 aa = &a[k3 + 1]; 108 for ( ll=0; ll<i__2; ll++ ) aa[ll] *= stmp; 109 kp1 = k + 1; 110 if (2 < kp1) continue; 111 ax = aa; 112 for (j = kp1; j <= 2; ++j) { 113 j3 = 2*j; 114 stmp = a[k + j3]; 115 a[k + j3] = 0.0; 116 ay = &a[j3 + 1]; 117 for ( ll=0; ll<k; ll++ ) { 118 ay[ll] += stmp*ax[ll]; 119 } 120 } 121 } 122 123 /* form inverse(u)*inverse(l) */ 124 125 /*for (kb = 1; kb <= 1; ++kb) {*/ 126 127 k = 1; 128 k3 = 2*k; 129 kp1 = k + 1; 130 aa = a + k3; 131 for (i = kp1; i <= 2; ++i) { 132 work_l[i-1] = aa[i]; 133 /* work[i] = aa[i]; Fix for -O3 error on Origin 2000 */ 134 aa[i] = 0.0; 135 } 136 for (j = kp1; j <= 2; ++j) { 137 stmp = work[j]; 138 ax = &a[2*j + 1]; 139 ay = &a[k3 + 1]; 140 ay[0] += stmp*ax[0]; 141 ay[1] += stmp*ax[1]; 142 } 143 l = ipvt[k]; 144 if (l != k) { 145 ax = &a[k3 + 1]; 146 ay = &a[2*l + 1]; 147 stmp = ax[0]; ax[0] = ay[0]; ay[0] = stmp; 148 stmp = ax[1]; ax[1] = ay[1]; ay[1] = stmp; 149 } 150 151 PetscFunctionReturn(0); 152 } 153 154 155