1be1d678aSKris Buschelman 2a97b7df5SSatish Balay /* 39fd41dd0SBarry Smith This routine was converted by f2c from Linpack source 49fd41dd0SBarry Smith linpack. this version dated 08/14/78 59fd41dd0SBarry Smith cleve moler, university of new mexico, argonne national lab. 659539b86SBarry Smith 759539b86SBarry Smith Does an LU factorization with partial pivoting of a dense 89d8b60e5SBarry Smith n by n matrix. 99d8b60e5SBarry Smith 109d8b60e5SBarry Smith Used by the sparse factorization routines in 11dd882469SBarry Smith src/mat/impls/baij/seq 1259539b86SBarry Smith 139fd41dd0SBarry Smith */ 14c6db04a5SJed Brown #include <petscsys.h> 159fd41dd0SBarry Smith 16d71ae5a4SJacob Faibussowitsch PETSC_INTERN PetscErrorCode PetscLINPACKgefa(MatScalar *a, PetscInt n, PetscInt *ipvt, PetscBool allowzeropivot, PetscBool *zeropivotdetected) 17d71ae5a4SJacob Faibussowitsch { 18690b6cddSBarry Smith PetscInt i__2, i__3, kp1, nm1, j, k, l, ll, kn, knp1, jn1; 193f1db9ecSBarry Smith MatScalar t, *ax, *ay, *aa; 20329f5518SBarry Smith MatReal tmp, max; 219fd41dd0SBarry Smith 223a40ed3dSBarry Smith PetscFunctionBegin; 235f8bbccaSHong Zhang if (zeropivotdetected) *zeropivotdetected = PETSC_FALSE; 245f8bbccaSHong Zhang 259fd41dd0SBarry Smith /* Parameter adjustments */ 269fd41dd0SBarry Smith --ipvt; 2739d66777SBarry Smith a -= n + 1; 289fd41dd0SBarry Smith 299fd41dd0SBarry Smith /* Function Body */ 309fd41dd0SBarry Smith nm1 = n - 1; 3139d66777SBarry Smith for (k = 1; k <= nm1; ++k) { 329fd41dd0SBarry Smith kp1 = k + 1; 3339d66777SBarry Smith kn = k * n; 3439d66777SBarry Smith knp1 = k * n + k; 359fd41dd0SBarry Smith 369fd41dd0SBarry Smith /* find l = pivot index */ 379fd41dd0SBarry Smith 389fd41dd0SBarry Smith i__2 = n - k + 1; 3939d66777SBarry Smith aa = &a[knp1]; 409fd41dd0SBarry Smith max = PetscAbsScalar(aa[0]); 419fd41dd0SBarry Smith l = 1; 429fd41dd0SBarry Smith for (ll = 1; ll < i__2; ll++) { 439fd41dd0SBarry Smith tmp = PetscAbsScalar(aa[ll]); 449371c9d4SSatish Balay if (tmp > max) { 459371c9d4SSatish Balay max = tmp; 469371c9d4SSatish Balay l = ll + 1; 479371c9d4SSatish Balay } 489fd41dd0SBarry Smith } 499fd41dd0SBarry Smith l += k - 1; 509fd41dd0SBarry Smith ipvt[k] = l; 519fd41dd0SBarry Smith 525f8bbccaSHong Zhang if (a[l + kn] == 0.0) { 535f8bbccaSHong Zhang if (allowzeropivot) { 549566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL, "Zero pivot, row %" PetscInt_FMT "\n", k - 1)); 5570d8d27cSBarry Smith if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE; 5698921bdaSJacob Faibussowitsch } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT, k - 1); 575f8bbccaSHong Zhang } 589fd41dd0SBarry Smith 599fd41dd0SBarry Smith /* interchange if necessary */ 609fd41dd0SBarry Smith if (l != k) { 6139d66777SBarry Smith t = a[l + kn]; 6239d66777SBarry Smith a[l + kn] = a[knp1]; 6339d66777SBarry Smith a[knp1] = t; 649fd41dd0SBarry Smith } 659fd41dd0SBarry Smith 669fd41dd0SBarry Smith /* compute multipliers */ 6739d66777SBarry Smith t = -1. / a[knp1]; 689fd41dd0SBarry Smith i__2 = n - k; 6939d66777SBarry Smith aa = &a[1 + knp1]; 7026fbe8dcSKarl Rupp for (ll = 0; ll < i__2; ll++) aa[ll] *= t; 719fd41dd0SBarry Smith 729fd41dd0SBarry Smith /* row elimination with column indexing */ 7339d66777SBarry Smith ax = aa; 749fd41dd0SBarry Smith for (j = kp1; j <= n; ++j) { 758d3e6ddaSBarry Smith jn1 = j * n; 768d3e6ddaSBarry Smith t = a[l + jn1]; 779fd41dd0SBarry Smith if (l != k) { 788d3e6ddaSBarry Smith a[l + jn1] = a[k + jn1]; 798d3e6ddaSBarry Smith a[k + jn1] = t; 809fd41dd0SBarry Smith } 819fd41dd0SBarry Smith 829fd41dd0SBarry Smith i__3 = n - k; 838d3e6ddaSBarry Smith ay = &a[1 + k + jn1]; 8426fbe8dcSKarl Rupp for (ll = 0; ll < i__3; ll++) ay[ll] += t * ax[ll]; 859fd41dd0SBarry Smith } 869fd41dd0SBarry Smith } 879fd41dd0SBarry Smith ipvt[n] = n; 885f8bbccaSHong Zhang if (a[n + n * n] == 0.0) { 895f8bbccaSHong Zhang if (allowzeropivot) { 909566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL, "Zero pivot, row %" PetscInt_FMT "\n", n - 1)); 9170d8d27cSBarry Smith if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE; 9298921bdaSJacob Faibussowitsch } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT, n - 1); 935f8bbccaSHong Zhang } 94*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 959fd41dd0SBarry Smith } 96