/* dgtsl.c --- LINPACK solver for general tridiagonal matrix. */ #include #include #include /* This is a translation of the FORTRAN routine from LINPACK. * * INPUT parameters: * n - order of tridiagonal matrix, * c - singly subscripted array with subdiagonal in * positions c(1) to c(n-1), * d - singly subscripted array with diagonal in * positions d(0) to d(n-1), * e - singly subscripted array with superdiagonal in * positions e(0) to e(n-2), * b - singly subscripted array with right hand side * of system, * * RETURN values: * The routine returns an integer which is 0 on successful * completion and contains the index+1 of the first zero * pivot found during the solution process otherwise. * b - solution vector, */ int dgtsl(double c[],double d[],double e[],double b[],int n) { double t; int nm1,k,kb,kp1,nm2,retval; /* Initialize return value. */ retval = 0; /* Use C array indexing rules. */ n--; c[0] = d[0]; nm1 = n-1; if (n > 0) { d[0] = e[0]; e[0] = 0.0; e[n] = 0.0; for (k = 0; k <= nm1; k++) { kp1 = k + 1; /* find the largest of the two rows */ if (fabs(c[kp1]) >= fabs(c[k])) { /* interchange rows */ t = c[kp1]; c[kp1] = c[k]; c[k] = t; t = d[kp1]; d[kp1] = d[k]; d[k] = t; t = e[kp1]; e[kp1] = e[k]; e[k] = t; t = b[kp1]; b[kp1] = b[k]; b[k] = t; } /* zero elements */ if (c[k] == 0.0) return k+1; /* .............exit */ t = -c[kp1] / c[k]; c[kp1] = d[kp1] + t * d[k]; d[kp1] = e[kp1] + t * e[k]; e[kp1] = 0.0; b[kp1] = b[kp1] + t * b[k]; } } if (c[n] == 0.0) return n+1; /* back solve */ nm2 = n - 2; b[n] = b[n] / c[n]; if (n != 0) { b[nm1] = (b[nm1] - d[nm1] * b[n]) / c[nm1]; if (nm2 >= 0) { for (kb = 0; kb <= nm2; kb++) { k = nm2 - kb; b[k] = (b[k] - d[k] * b[k+1] - e[k] * b[k+2]) / c[k]; } /* endif */ } /* endif */ } return retval; }