1 2 /* 3 Inverts 2 by 2 matrix using partial pivoting. 4 5 Used by the sparse factorization routines in 6 src/mat/impls/baij/seq 7 8 9 This is a combination of the Linpack routines 10 dgefa() and dgedi() specialized for a size of 2. 11 12 */ 13 #include <petscsys.h> 14 15 #undef __FUNCT__ 16 #define __FUNCT__ "PetscKernel_A_gets_inverse_A_2" 17 PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_2(MatScalar *a,PetscReal shift,PetscBool allowzeropivot,PetscBool *zeropivotdetected) 18 { 19 PetscInt i__2,i__3,kp1,j,k,l,ll,i,ipvt[2],k3; 20 PetscInt k4,j3; 21 MatScalar *aa,*ax,*ay,work[4],stmp; 22 MatReal tmp,max; 23 24 /* gaussian elimination with partial pivoting */ 25 26 PetscFunctionBegin; 27 shift = .25*shift*(1.e-12 + PetscAbsScalar(a[0]) + PetscAbsScalar(a[3])); 28 /* Parameter adjustments */ 29 a -= 3; 30 31 /*for (k = 1; k <= 1; ++k) {*/ 32 k = 1; 33 kp1 = k + 1; 34 k3 = 2*k; 35 k4 = k3 + k; 36 /* find l = pivot index */ 37 38 i__2 = 3 - k; 39 aa = &a[k4]; 40 max = PetscAbsScalar(aa[0]); 41 l = 1; 42 for (ll=1; ll<i__2; ll++) { 43 tmp = PetscAbsScalar(aa[ll]); 44 if (tmp > max) { max = tmp; l = ll+1;} 45 } 46 l += k - 1; 47 ipvt[k-1] = l; 48 49 if (a[l + k3] == 0.0) { 50 if (shift == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",k-1); 51 else { 52 a[l + k3] = shift; 53 } 54 } 55 56 /* interchange if necessary */ 57 58 if (l != k) { 59 stmp = a[l + k3]; 60 a[l + k3] = a[k4]; 61 a[k4] = stmp; 62 } 63 64 /* compute multipliers */ 65 66 stmp = -1. / a[k4]; 67 i__2 = 2 - k; 68 aa = &a[1 + k4]; 69 for (ll=0; ll<i__2; ll++) aa[ll] *= stmp; 70 71 /* row elimination with column indexing */ 72 73 ax = &a[k4+1]; 74 for (j = kp1; j <= 2; ++j) { 75 j3 = 2*j; 76 stmp = a[l + j3]; 77 if (l != k) { 78 a[l + j3] = a[k + j3]; 79 a[k + j3] = stmp; 80 } 81 82 i__3 = 2 - k; 83 ay = &a[1+k+j3]; 84 for (ll=0; ll<i__3; ll++) ay[ll] += stmp*ax[ll]; 85 } 86 87 ipvt[1] = 2; 88 if (a[6] == 0.0) { 89 PetscErrorCode ierr; 90 if (allowzeropivot) { 91 ierr = PetscInfo1(NULL,"Zero pivot, row %D\n",1);CHKERRQ(ierr); 92 *zeropivotdetected = PETSC_TRUE; 93 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",1); 94 } 95 96 /* 97 Now form the inverse 98 */ 99 100 /* compute inverse(u) */ 101 102 for (k = 1; k <= 2; ++k) { 103 k3 = 2*k; 104 k4 = k3 + k; 105 a[k4] = 1.0 / a[k4]; 106 stmp = -a[k4]; 107 i__2 = k - 1; 108 aa = &a[k3 + 1]; 109 for (ll=0; ll<i__2; ll++) aa[ll] *= stmp; 110 kp1 = k + 1; 111 if (2 < kp1) continue; 112 ax = aa; 113 for (j = kp1; j <= 2; ++j) { 114 j3 = 2*j; 115 stmp = a[k + j3]; 116 a[k + j3] = 0.0; 117 ay = &a[j3 + 1]; 118 for (ll=0; ll<k; ll++) ay[ll] += stmp*ax[ll]; 119 } 120 } 121 122 /* form inverse(u)*inverse(l) */ 123 k = 1; 124 k3 = 2*k; 125 kp1 = k + 1; 126 aa = a + k3; 127 for (i = kp1; i <= 2; ++i) { 128 work[i-1] = aa[i]; 129 aa[i] = 0.0; 130 } 131 for (j = kp1; j <= 2; ++j) { 132 stmp = work[j-1]; 133 ax = &a[2*j + 1]; 134 ay = &a[k3 + 1]; 135 ay[0] += stmp*ax[0]; 136 ay[1] += stmp*ax[1]; 137 } 138 l = ipvt[k-1]; 139 if (l != k) { 140 ax = &a[k3 + 1]; 141 ay = &a[2*l + 1]; 142 stmp = ax[0]; ax[0] = ay[0]; ay[0] = stmp; 143 stmp = ax[1]; ax[1] = ay[1]; ay[1] = stmp; 144 } 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNCT__ 149 #define __FUNCT__ "PetscKernel_A_gets_inverse_A_9" 150 PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_9(MatScalar *a,PetscReal shift,PetscBool allowzeropivot,PetscBool *zeropivotdetected) 151 { 152 PetscInt i__2,i__3,kp1,j,k,l,ll,i,ipvt[9],kb,k3; 153 PetscInt k4,j3; 154 MatScalar *aa,*ax,*ay,work[81],stmp; 155 MatReal tmp,max; 156 157 /* gaussian elimination with partial pivoting */ 158 159 PetscFunctionBegin; 160 /* Parameter adjustments */ 161 a -= 10; 162 163 for (k = 1; k <= 8; ++k) { 164 kp1 = k + 1; 165 k3 = 9*k; 166 k4 = k3 + k; 167 /* find l = pivot index */ 168 169 i__2 = 10 - k; 170 aa = &a[k4]; 171 max = PetscAbsScalar(aa[0]); 172 l = 1; 173 for (ll=1; ll<i__2; ll++) { 174 tmp = PetscAbsScalar(aa[ll]); 175 if (tmp > max) { max = tmp; l = ll+1;} 176 } 177 l += k - 1; 178 ipvt[k-1] = l; 179 180 if (a[l + k3] == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",k-1); 181 182 /* interchange if necessary */ 183 184 if (l != k) { 185 stmp = a[l + k3]; 186 a[l + k3] = a[k4]; 187 a[k4] = stmp; 188 } 189 190 /* compute multipliers */ 191 192 stmp = -1. / a[k4]; 193 i__2 = 9 - k; 194 aa = &a[1 + k4]; 195 for (ll=0; ll<i__2; ll++) aa[ll] *= stmp; 196 197 /* row elimination with column indexing */ 198 199 ax = &a[k4+1]; 200 for (j = kp1; j <= 9; ++j) { 201 j3 = 9*j; 202 stmp = a[l + j3]; 203 if (l != k) { 204 a[l + j3] = a[k + j3]; 205 a[k + j3] = stmp; 206 } 207 208 i__3 = 9 - k; 209 ay = &a[1+k+j3]; 210 for (ll=0; ll<i__3; ll++) ay[ll] += stmp*ax[ll]; 211 } 212 } 213 ipvt[8] = 9; 214 if (a[90] == 0.0) { 215 PetscErrorCode ierr; 216 if (allowzeropivot) { 217 ierr = PetscInfo1(NULL,"Zero pivot, row %D\n",6);CHKERRQ(ierr); 218 *zeropivotdetected = PETSC_TRUE; 219 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",6); 220 } 221 222 /* 223 Now form the inverse 224 */ 225 226 /* compute inverse(u) */ 227 228 for (k = 1; k <= 9; ++k) { 229 k3 = 9*k; 230 k4 = k3 + k; 231 a[k4] = 1.0 / a[k4]; 232 stmp = -a[k4]; 233 i__2 = k - 1; 234 aa = &a[k3 + 1]; 235 for (ll=0; ll<i__2; ll++) aa[ll] *= stmp; 236 kp1 = k + 1; 237 if (9 < kp1) continue; 238 ax = aa; 239 for (j = kp1; j <= 9; ++j) { 240 j3 = 9*j; 241 stmp = a[k + j3]; 242 a[k + j3] = 0.0; 243 ay = &a[j3 + 1]; 244 for (ll=0; ll<k; ll++) ay[ll] += stmp*ax[ll]; 245 } 246 } 247 248 /* form inverse(u)*inverse(l) */ 249 250 for (kb = 1; kb <= 8; ++kb) { 251 k = 9 - kb; 252 k3 = 9*k; 253 kp1 = k + 1; 254 aa = a + k3; 255 for (i = kp1; i <= 9; ++i) { 256 work[i-1] = aa[i]; 257 aa[i] = 0.0; 258 } 259 for (j = kp1; j <= 9; ++j) { 260 stmp = work[j-1]; 261 ax = &a[9*j + 1]; 262 ay = &a[k3 + 1]; 263 ay[0] += stmp*ax[0]; 264 ay[1] += stmp*ax[1]; 265 ay[2] += stmp*ax[2]; 266 ay[3] += stmp*ax[3]; 267 ay[4] += stmp*ax[4]; 268 ay[5] += stmp*ax[5]; 269 ay[6] += stmp*ax[6]; 270 ay[7] += stmp*ax[7]; 271 ay[8] += stmp*ax[8]; 272 } 273 l = ipvt[k-1]; 274 if (l != k) { 275 ax = &a[k3 + 1]; 276 ay = &a[9*l + 1]; 277 stmp = ax[0]; ax[0] = ay[0]; ay[0] = stmp; 278 stmp = ax[1]; ax[1] = ay[1]; ay[1] = stmp; 279 stmp = ax[2]; ax[2] = ay[2]; ay[2] = stmp; 280 stmp = ax[3]; ax[3] = ay[3]; ay[3] = stmp; 281 stmp = ax[4]; ax[4] = ay[4]; ay[4] = stmp; 282 stmp = ax[5]; ax[5] = ay[5]; ay[5] = stmp; 283 stmp = ax[6]; ax[6] = ay[6]; ay[6] = stmp; 284 stmp = ax[7]; ax[7] = ay[7]; ay[7] = stmp; 285 stmp = ax[8]; ax[8] = ay[8]; ay[8] = stmp; 286 } 287 } 288 PetscFunctionReturn(0); 289 } 290 291 /* 292 Inverts 15 by 15 matrix using partial pivoting. 293 294 Used by the sparse factorization routines in 295 src/mat/impls/baij/seq 296 297 This is a combination of the Linpack routines 298 dgefa() and dgedi() specialized for a size of 15. 299 300 */ 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "PetscKernel_A_gets_inverse_A_15" 304 PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_15(MatScalar *a,PetscInt *ipvt,MatScalar *work,PetscReal shift) 305 { 306 PetscInt i__2,i__3,kp1,j,k,l,ll,i,kb,k3; 307 PetscInt k4,j3; 308 MatScalar *aa,*ax,*ay,stmp; 309 MatReal tmp,max; 310 311 /* gaussian elimination with partial pivoting */ 312 313 PetscFunctionBegin; 314 /* Parameter adjustments */ 315 a -= 16; 316 317 for (k = 1; k <= 14; ++k) { 318 kp1 = k + 1; 319 k3 = 15*k; 320 k4 = k3 + k; 321 /* find l = pivot index */ 322 323 i__2 = 16 - k; 324 aa = &a[k4]; 325 max = PetscAbsScalar(aa[0]); 326 l = 1; 327 for (ll=1; ll<i__2; ll++) { 328 tmp = PetscAbsScalar(aa[ll]); 329 if (tmp > max) { max = tmp; l = ll+1;} 330 } 331 l += k - 1; 332 ipvt[k-1] = l; 333 334 if (a[l + k3] == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",k-1); 335 336 /* interchange if necessary */ 337 338 if (l != k) { 339 stmp = a[l + k3]; 340 a[l + k3] = a[k4]; 341 a[k4] = stmp; 342 } 343 344 /* compute multipliers */ 345 346 stmp = -1. / a[k4]; 347 i__2 = 15 - k; 348 aa = &a[1 + k4]; 349 for (ll=0; ll<i__2; ll++) aa[ll] *= stmp; 350 351 /* row elimination with column indexing */ 352 353 ax = &a[k4+1]; 354 for (j = kp1; j <= 15; ++j) { 355 j3 = 15*j; 356 stmp = a[l + j3]; 357 if (l != k) { 358 a[l + j3] = a[k + j3]; 359 a[k + j3] = stmp; 360 } 361 362 i__3 = 15 - k; 363 ay = &a[1+k+j3]; 364 for (ll=0; ll<i__3; ll++) ay[ll] += stmp*ax[ll]; 365 } 366 } 367 ipvt[14] = 15; 368 if (a[240] == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",6); 369 370 /* 371 Now form the inverse 372 */ 373 374 /* compute inverse(u) */ 375 376 for (k = 1; k <= 15; ++k) { 377 k3 = 15*k; 378 k4 = k3 + k; 379 a[k4] = 1.0 / a[k4]; 380 stmp = -a[k4]; 381 i__2 = k - 1; 382 aa = &a[k3 + 1]; 383 for (ll=0; ll<i__2; ll++) aa[ll] *= stmp; 384 kp1 = k + 1; 385 if (15 < kp1) continue; 386 ax = aa; 387 for (j = kp1; j <= 15; ++j) { 388 j3 = 15*j; 389 stmp = a[k + j3]; 390 a[k + j3] = 0.0; 391 ay = &a[j3 + 1]; 392 for (ll=0; ll<k; ll++) ay[ll] += stmp*ax[ll]; 393 } 394 } 395 396 /* form inverse(u)*inverse(l) */ 397 398 for (kb = 1; kb <= 14; ++kb) { 399 k = 15 - kb; 400 k3 = 15*k; 401 kp1 = k + 1; 402 aa = a + k3; 403 for (i = kp1; i <= 15; ++i) { 404 work[i-1] = aa[i]; 405 aa[i] = 0.0; 406 } 407 for (j = kp1; j <= 15; ++j) { 408 stmp = work[j-1]; 409 ax = &a[15*j + 1]; 410 ay = &a[k3 + 1]; 411 ay[0] += stmp*ax[0]; 412 ay[1] += stmp*ax[1]; 413 ay[2] += stmp*ax[2]; 414 ay[3] += stmp*ax[3]; 415 ay[4] += stmp*ax[4]; 416 ay[5] += stmp*ax[5]; 417 ay[6] += stmp*ax[6]; 418 ay[7] += stmp*ax[7]; 419 ay[8] += stmp*ax[8]; 420 ay[9] += stmp*ax[9]; 421 ay[10] += stmp*ax[10]; 422 ay[11] += stmp*ax[11]; 423 ay[12] += stmp*ax[12]; 424 ay[13] += stmp*ax[13]; 425 ay[14] += stmp*ax[14]; 426 } 427 l = ipvt[k-1]; 428 if (l != k) { 429 ax = &a[k3 + 1]; 430 ay = &a[15*l + 1]; 431 stmp = ax[0]; ax[0] = ay[0]; ay[0] = stmp; 432 stmp = ax[1]; ax[1] = ay[1]; ay[1] = stmp; 433 stmp = ax[2]; ax[2] = ay[2]; ay[2] = stmp; 434 stmp = ax[3]; ax[3] = ay[3]; ay[3] = stmp; 435 stmp = ax[4]; ax[4] = ay[4]; ay[4] = stmp; 436 stmp = ax[5]; ax[5] = ay[5]; ay[5] = stmp; 437 stmp = ax[6]; ax[6] = ay[6]; ay[6] = stmp; 438 stmp = ax[7]; ax[7] = ay[7]; ay[7] = stmp; 439 stmp = ax[8]; ax[8] = ay[8]; ay[8] = stmp; 440 stmp = ax[9]; ax[9] = ay[9]; ay[9] = stmp; 441 stmp = ax[10]; ax[10] = ay[10]; ay[10] = stmp; 442 stmp = ax[11]; ax[11] = ay[11]; ay[11] = stmp; 443 stmp = ax[12]; ax[12] = ay[12]; ay[12] = stmp; 444 stmp = ax[13]; ax[13] = ay[13]; ay[13] = stmp; 445 stmp = ax[14]; ax[14] = ay[14]; ay[14] = stmp; 446 } 447 } 448 PetscFunctionReturn(0); 449 } 450