Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/lu.cpp
Go to the documentation of this file.
00001 /******************************************************
00002 *    LU decomposition routines used by test_lu.cpp    *
00003 *    with dynamic allocations                         *
00004 *                                                     *
00005 *                  C++ version by J-P Moreau, Paris   *
00006 * --------------------------------------------------- *
00007 * Reference:                                          *
00008 *                                                     *
00009 * "Numerical Recipes by W.H. Press, B. P. Flannery,   *
00010 *  S.A. Teukolsky and W.T. Vetterling, Cambridge      *
00011 *  University Press, 1986".                           *
00012 * --------------------------------------------------- *
00013 * Uses: basis_r.cpp and vmblock.cpp                   *
00014 ******************************************************/
00015 
00016 #include "lu.h"
00017 
00018 /**************************************************************
00019 * Given an N x N matrix A, this routine replaces it by the LU *
00020 * decomposition of a rowwise permutation of itself. A and N   *
00021 * are input. INDX is an output vector which records the row   *
00022 * permutation effected by the partial pivoting; D is output   *
00023 * as -1 or 1, depending on whether the number of row inter-   *
00024 * changes was even or odd, respectively. This routine is used *
00025 * in combination with LUBKSB to solve linear equations or to  *
00026 * invert a matrix. Return code is 1, if matrix is singular.   *
00027 **************************************************************/
00028 int LUDCMP(REAL **A, int n, int *INDX, int *d)  {
00029 
00030   REAL AMAX,DUM, SUM;
00031   int  I,IMAX,J,K;
00032   IMAX = 0;
00033 
00034   REAL *VV;
00035   void *vmblock = NULL;
00036 
00037   vmblock = vminit();
00038   VV      = (REAL *) vmalloc(vmblock, VEKTOR,  NMAX, 0);
00039 
00040   if (! vmcomplete(vmblock))
00041   {
00042     // LogError ("No Memory", 0, __FILE__, __LINE__);
00043     return -1;
00044   }
00045 
00046   *d=1;
00047 
00048   for  (I=1; I<n+1; I++)  {
00049     AMAX=0.0;
00050     for (J=1; J<n+1; J++)
00051       if (ABS(A[I][J]) > AMAX)  AMAX=ABS(A[I][J]);
00052 
00053     if(AMAX < TINY)
00054       return 1; //singular Matrix!
00055     VV[I] = 1.0 / AMAX;
00056   } // i loop
00057 
00058   for (J=1; J<n+1; J++)  {
00059 
00060     for (I=1; I<J; I++)  {
00061       SUM = A[I][J];
00062       for (K=1; K<I; K++)
00063         SUM = SUM - A[I][K]*A[K][J];
00064       A[I][J] = SUM;
00065     } // i loop
00066     AMAX = 0.0;
00067 
00068     for (I=J; I<n+1; I++)  {
00069       SUM = A[I][J];
00070       for  (K=1; K<J; K++)
00071         SUM = SUM - A[I][K]*A[K][J];
00072       A[I][J] = SUM;
00073       DUM = VV[I]*ABS(SUM);
00074       if (DUM >= AMAX) {
00075         IMAX = I;
00076         AMAX = DUM;
00077       }
00078     } // i loop
00079 
00080     if (J != IMAX)  {
00081       for (K=1; K<n+1; K++)  {
00082         DUM = A[IMAX][K];
00083         A[IMAX][K] = A[J][K];
00084         A[J][K] = DUM;
00085       } // k loop
00086       *d = -*d;
00087       VV[IMAX] = VV[J];
00088     }
00089 
00090     INDX[J] = IMAX;
00091 
00092     if (ABS(A[J][J]) < TINY)   A[J][J] = TINY;
00093 
00094     if (J != n)  {
00095       DUM = 1.0 / A[J][J];
00096       for (I=J+1; I<n+1; I++)
00097         A[I][J] *= DUM;
00098     }
00099   } // j loop
00100 
00101   free(vmblock);
00102   return 0;
00103 
00104 } // subroutine LUDCMP
00105 
00106 
00107 /*****************************************************************
00108 * Solves the set of N linear equations A . X = B.  Here A is     *
00109 * input, not as the matrix A but rather as its LU decomposition, *
00110 * determined by the routine LUDCMP. INDX is input as the permuta-*
00111 * tion vector returned by LUDCMP. B is input as the right-hand   *
00112 * side vector B, and returns with the solution vector X. A, N and*
00113 * INDX are not modified by this routine and can be used for suc- *
00114 * cessive calls with different right-hand sides. This routine is *
00115 * also efficient for plain matrix inversion.                     *
00116 *****************************************************************/
00117 void LUBKSB(REAL **A, int n, int *INDX, REAL *B)  {
00118   REAL SUM;
00119   int  I,II,J,LL;
00120 
00121   II = 0;
00122 
00123   for (I=1; I<n+1; I++)  {
00124     LL = INDX[I];
00125     SUM = B[LL];
00126     B[LL] = B[I];
00127     if (II != 0)
00128       for (J=II; J<I; J++)
00129         SUM = SUM - A[I][J]*B[J];
00130     else if (SUM != 0.0)  II = I;
00131     B[I] = SUM;
00132   } // i loop
00133 
00134   for (I=n; I>0; I--)  {
00135     SUM = B[I];
00136     if (I < n)  {
00137       for (J=I+1; J<n+1; J++)
00138         SUM = SUM - A[I][J]*B[J];
00139     }
00140     B[I] = SUM / A[I][I];
00141   } // i loop
00142 
00143 
00144 } // LUBKSB
00145 
00146 // end of file lu.cpp
00147 
00148