#include #include //================================================================C int GLU3(double *A, double *B, int N, int ND, double EPS, int *IP) //================================================================C // Real-Dense LU Decomposition by Gauss Elimination C // A ---> LU Decomposition with Partial Pivoting C // Type-1 Partial Pivoting (Changing All rows) C //----------------------------------------------------------------C // A(ND*N) R*8, I/O, Coefficient Matrix C // B(N) R*8, I/O, Righ-hand Vector C // N I*4, In, Matrix Size of A C // ND I*4, In, Array Size of A ( ND >= N ) C // EPS R*8, In, Value for Singularity Check C // IP(N) I*4, Out, Pivot Number C // return I*4, Out, 0 : Normal Execution C // 1 : Singular Stop C // 2 ; Parameter Error C //----------------------------------------------------------------C // Written by Yasunori Ushiro , 2007/05/30 C // ( Tokyo Polytechnic University ) C // 後 保範(東京工芸大学) C // Ver.1 : No tuning Version ; Address is equal to Fortran C //================================================================C { int i, j, k, KPIV, IER, KND, JND ; double PIV, AKJ, DPIV, BW ; //----- Gauss Elimination ----------------------------- // Check Parameter IER = 0 ; if(ND < N) { IER = 2 ; goto M100 ; } // Gauss Elimination for (k=0; k= PIV) { KPIV = i ; PIV = fabs(A[i+KND]) ; } } // Check Singularity if(PIV < EPS) { IER = 1 ; goto M100 ; } IP[k] = KPIV ; // Change A(k,k) <--> A(KPIV,k) if(KPIV != k) { for (j=0; j=0; j--) { B[j] = B[j]*A[j+j*ND] ; for (i=0; i<=j-1; i++) { B[i] = B[i] - A[i+j*ND]*B[j] ; } } M100: return (IER) ; }