/************************************************************* * Solve banded linear system AX = B By LU decomposition * * ---------------------------------------------------------- * * SAMPLE RUN: * * (Solve banded linear system AX = B, where: * * 1 10 0 0 0 1 * * 9 2 20 0 0 1 * * A = 0 19 3 30 0 B = 1 * * 0 0 29 4 40 1 * * 0 0 0 39 5 1 ) * * * * Note: A is given as follows: * * 0 0 0 0 0 * * 0 10 20 30 40 * * 1 2 3 4 5 * * 9 19 29 39 0 * * 0 0 0 0 0 * * * * SOLUTION: * * * * 0.40664649 * * 0.05933535 * * -0.13892446 * * 0.00964672 * * 0.12475556 * * * * Error code: 0 * * * * ---------------------------------------------------------- * * Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * * [BIBLI 18]. * * * * C++ Release 1.0 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *************************************************************/ #include #include #define NMAX 25 typedef double MAT[NMAX][NMAX]; MAT A; //banded matrix double B[NMAX]; //second member vector double X[NMAX]; //solution vector int IPVT[NMAX]; //auxiliary integer vector int IND; //error code int I,J,N,MU,ML; int Max(int a, int b) { if (a>=b) return a; else return b; } int Min(int a, int b) { if (a<=b) return a; else return b; } int IAMAX(double *A, int N) { double T; int I; T=0.0; for (I=1; I<=N; I++) if (fabs(A[I]) > T) { T=fabs(A[I]); return I; } return 0; } void SCALE(int N, double T, double *A) { int I; for (I=1; I<=N; I++) A[I]=T*A[I]; } void DAXPY(int N, double A, double *X, double *Y) { int I; for (I=1; I<=N; I++) Y[I] = Y[I] + A*X[I]; } void NSBFAC(MAT B, int N, int ML, int MU, int *IPVT, int *IND) { /*------------------------------------------------------------------ ! LU factorization of a band matrix (non symmetric) with partial ! pivoting. ! INPUTS: ! B : banded matrix. The correspondance between full matrix ! A(i,j) and band matrix B(k,l) is given by following sequence: ! m=ml+mu+1; ! for (j=1; j<=n; j++) { ! i1=Max(1,j-mu); ! i2=Min(n,j+ml); ! for (i=i1; i<=i2; i++) { ! k=i-j+m; ! B[k][j] = A[i][j]; ! } ! } ! N : size of B ! ML : number of lower diagonals ! MU : number of upper diagonals ! OUTPUTS: ! B : banded matrix storing the LU elements, the ML first lines ! of which are used for pivoting. ! IPVT: integer vector of size N storing the pivoting indices. ! IND : flag = 0, B is non singular ! = k, B may be singular !------------------------------------------------------------------*/ // Labels: e10, e20; double T; int I0,I,J,J0,J1,JU,JZ,K,KP1,L,LM,M,MM,NM1; double TMP[NMAX], TMP1[NMAX]; M=ML+MU+1; *IND=0; J0=MU+2; J1=Min(N,M)-1; if (J1 >= J0) for (JZ=J0; JZ<=J1; JZ++) { I0=M+1-JZ; for (I=I0; I<=ML; I++) B[I][JZ]=0.0; } JZ=J1; JU=0; NM1=N-1; if (NM1 >= 1) for (K=1; K<=NM1; K++) { KP1=K+1; JZ=JZ+1; if (JZ <= N && ML >= 1) for (I=1; I<=ML; I++) B[I][JZ]=0.0; LM=Min(ML,N-K); for (I=1; I<=LM+1; I++) TMP[I]=B[I+M-1][K]; L=IAMAX(TMP,LM+1)+M-1; IPVT[K]=L+K-M; if (B[L][K] == 0.0) goto e10; if (L != M) { T=B[L][K]; B[L][K]=B[M][K]; B[M][K]=T; } T=-1.0/B[M][K]; for (I=1; I<=LM; I++) TMP[I]=B[I+M][K]; SCALE(LM,T,TMP); for (I=1; I<=LM; I++) B[I+M][K]=TMP[I]; JU=Min(Max(JU,MU+IPVT[K]),N); MM=M; if (JU >= KP1) for (J=KP1; J<=JU; J++) { L--; MM--; T=B[L][J]; if (L != MM) { B[L][J]=B[MM][J]; B[MM][J]=T; } for (I=1; I<=LM; I++) { TMP[I]=B[I+M][K]; TMP1[I]=B[I+MM][J]; } DAXPY(LM,T,TMP,TMP1); for (I=1; I<=LM; I++) B[I+MM][J]=TMP1[I]; } goto e20; e10: *IND=K; e20: ;} IPVT[N]=N; if (B[M][N] == 0.0) *IND=N; } void NSBSLV(MAT A,int N,int ML,int MU, int *IPVT, double *B, double *X) { /*------------------------------------------------------------------- ! Solve banded linear system Ax := b ! INPUTS: ! A : banded matrix as output of LU factorization by NSBFAC ! (see storing mode in NSBFAC subroutine). ! LDA : 1st dimension of A in calling program (lda >= 2ml+mu+1) ! N : order of A --------------- ! ML : number of lower diagonals ! MU : number of upper diagonals ! IPVT: integer vector of size N storing the pivoting indices ! as output of NSBFAC. ! B : second member vector ! OUTPUT: ! X : solution vector !-------------------------------------------------------------------*/ double T; int I,K,L,LA,LB,LM,M,NM1; double TMP[NMAX],TMP1[NMAX]; for (I=1; I<=N; I++) X[I]=B[I]; M=ML+MU+1; NM1=N-1; // solve L*y = b if (ML != 0 && NM1 >= 1) for (K=1; K<=NM1; K++) { LM=Min(ML,N-K); L=IPVT[K]; T=X[L]; if (L != K) { X[L]=X[K]; X[K]=T; } for (I=1; I<=LM; I++) { TMP[I]=A[I+M][K]; TMP1[I]=X[I+K]; } DAXPY(LM,T,TMP,TMP1); for (I=1; I<=LM; I++) X[I+K]=TMP1[I]; } // solve U*y = x for (K=N; K>0; K--) { X[K]=X[K]/A[M][K]; LM=Min(K,M)-1; LA=M-LM; LB=K-LM; T=-X[K]; for (I=1; I<=LM; I++) { TMP[I]=A[I+LA-1][K]; TMP1[I]=X[I+LB-1]; } DAXPY(LM,T,TMP,TMP1); for (I=1; I<=LM; I++) X[I+LB-1]=TMP1[I]; } } void main() { N=5; //size of system for (I=1; I<=N; I++) for (J=1; J<=N; J++) A[I][J] = 0.0; MU=1; //one upper subdiagonal} ML=1; //one lower subdiagonal} // Define banded A matrix A[2][2]=10.0; A[2][3]=20.0; A[2][4]=30.0; A[2][5]=40.0; A[3][1]= 1.0; A[3][2]= 2.0; A[3][3]= 3.0; A[3][4]= 4.0; A[3][5]=5.0; A[4][1]= 9.0; A[4][2]=19.0; A[4][3]=29.0; A[4][4]=39.0; // Example #1 (B(1)..B(5) = ONE) for (I=1; I<=N; I++) B[I] = 1.0; // Example #2 (results are X1..X5 = ONE // B[1]=11.0; B[2]=31.0; B[3]=52.0; B[4]=73.0; B[5]=44.0; // call LU factorization NSBFAC(A,N,ML,MU,IPVT,&IND); // call appropriate solver NSBSLV(A,N,ML,MU,IPVT,B,X); // print solution and error code (must be zero) printf("\n SOLUTION:\n\n"); for (I=1; I<=N; I++) printf("%12.8f\n", X[I]); printf("\n Error code: %d\n\n", IND); } // end of file nsbslv.cpp