1*9fd41dd0SBarry Smith 2*9fd41dd0SBarry Smith /* 3*9fd41dd0SBarry Smith This routine was converted by f2c from Linpack source 4*9fd41dd0SBarry Smith linpack. this version dated 08/14/78 5*9fd41dd0SBarry Smith cleve moler, university of new mexico, argonne national lab. 6*9fd41dd0SBarry Smith */ 7*9fd41dd0SBarry Smith #include "petsc.h" 8*9fd41dd0SBarry Smith 9*9fd41dd0SBarry Smith int Linpack_DGEFA(Scalar *a, int n, int *ipvt) 10*9fd41dd0SBarry Smith { 11*9fd41dd0SBarry Smith int a_offset, i__1, i__2, i__3, kp1, nm1, j, k, l,ll; 12*9fd41dd0SBarry Smith Scalar t,*aa,tmp,max,*ax,*ay; 13*9fd41dd0SBarry Smith 14*9fd41dd0SBarry Smith /* gaussian elimination with partial pivoting */ 15*9fd41dd0SBarry Smith 16*9fd41dd0SBarry Smith /* Parameter adjustments */ 17*9fd41dd0SBarry Smith --ipvt; 18*9fd41dd0SBarry Smith a_offset = n + 1; 19*9fd41dd0SBarry Smith a -= a_offset; 20*9fd41dd0SBarry Smith 21*9fd41dd0SBarry Smith /* Function Body */ 22*9fd41dd0SBarry Smith nm1 = n - 1; 23*9fd41dd0SBarry Smith if (nm1 < 1) { 24*9fd41dd0SBarry Smith goto L70; 25*9fd41dd0SBarry Smith } 26*9fd41dd0SBarry Smith i__1 = nm1; 27*9fd41dd0SBarry Smith for (k = 1; k <= i__1; ++k) { 28*9fd41dd0SBarry Smith kp1 = k + 1; 29*9fd41dd0SBarry Smith 30*9fd41dd0SBarry Smith /* find l = pivot index */ 31*9fd41dd0SBarry Smith 32*9fd41dd0SBarry Smith i__2 = n - k + 1; 33*9fd41dd0SBarry Smith /* l = idamax_(&i__2, &a[k + k * n], &c__1) + k - 1; */ 34*9fd41dd0SBarry Smith aa = &a[k + k * n]; 35*9fd41dd0SBarry Smith max = PetscAbsScalar(aa[0]); 36*9fd41dd0SBarry Smith l = 1; 37*9fd41dd0SBarry Smith for ( ll=1; ll<i__2; ll++ ) { 38*9fd41dd0SBarry Smith tmp = PetscAbsScalar(aa[ll]); 39*9fd41dd0SBarry Smith if (tmp > max) { max = tmp; l = ll+1;} 40*9fd41dd0SBarry Smith } 41*9fd41dd0SBarry Smith l += k - 1; 42*9fd41dd0SBarry Smith ipvt[k] = l; 43*9fd41dd0SBarry Smith 44*9fd41dd0SBarry Smith /* zero pivot implies this column already triangularized */ 45*9fd41dd0SBarry Smith 46*9fd41dd0SBarry Smith if (a[l + k * n] == 0.) { 47*9fd41dd0SBarry Smith SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 48*9fd41dd0SBarry Smith } 49*9fd41dd0SBarry Smith 50*9fd41dd0SBarry Smith /* interchange if necessary */ 51*9fd41dd0SBarry Smith 52*9fd41dd0SBarry Smith if (l != k) { 53*9fd41dd0SBarry Smith t = a[l + k * n]; 54*9fd41dd0SBarry Smith a[l + k * n] = a[k + k * n]; 55*9fd41dd0SBarry Smith a[k + k * n] = t; 56*9fd41dd0SBarry Smith } 57*9fd41dd0SBarry Smith 58*9fd41dd0SBarry Smith /* compute multipliers */ 59*9fd41dd0SBarry Smith 60*9fd41dd0SBarry Smith t = -1. / a[k + k * n]; 61*9fd41dd0SBarry Smith i__2 = n - k; 62*9fd41dd0SBarry Smith /* dscal_(&i__2, &t, &a[k + 1 + k * n], &c__1); */ 63*9fd41dd0SBarry Smith aa = &a[k + 1 + k * n]; 64*9fd41dd0SBarry Smith for ( ll=0; ll<i__2; ll++ ) { 65*9fd41dd0SBarry Smith aa[ll] *= t; 66*9fd41dd0SBarry Smith } 67*9fd41dd0SBarry Smith 68*9fd41dd0SBarry Smith 69*9fd41dd0SBarry Smith /* row elimination with column indexing */ 70*9fd41dd0SBarry Smith 71*9fd41dd0SBarry Smith ax = &a[k+1+k*n]; 72*9fd41dd0SBarry Smith for (j = kp1; j <= n; ++j) { 73*9fd41dd0SBarry Smith t = a[l + j * n]; 74*9fd41dd0SBarry Smith if (l != k) { 75*9fd41dd0SBarry Smith a[l + j * n] = a[k + j * n]; 76*9fd41dd0SBarry Smith a[k + j * n] = t; 77*9fd41dd0SBarry Smith } 78*9fd41dd0SBarry Smith 79*9fd41dd0SBarry Smith i__3 = n - k; 80*9fd41dd0SBarry Smith /* daxpy_(&i__3, &t, &a[k + 1 + k * n], &c__1, &a[k + 1 + j * 81*9fd41dd0SBarry Smith n], &c__1); */ 82*9fd41dd0SBarry Smith ay = &a[k+1+j*n]; 83*9fd41dd0SBarry Smith for ( ll=0; ll<i__3; ll++ ) { 84*9fd41dd0SBarry Smith ay[ll] += t*ax[ll]; 85*9fd41dd0SBarry Smith } 86*9fd41dd0SBarry Smith } 87*9fd41dd0SBarry Smith } 88*9fd41dd0SBarry Smith L70: 89*9fd41dd0SBarry Smith ipvt[n] = n; 90*9fd41dd0SBarry Smith if (a[n + n * n] == 0.) { 91*9fd41dd0SBarry Smith SETERRQ(n,"Linpack_DGEFA:Zero pivot,final row"); 92*9fd41dd0SBarry Smith } 93*9fd41dd0SBarry Smith return 0; 94*9fd41dd0SBarry Smith } 95*9fd41dd0SBarry Smith 96