lu.c

Go to the documentation of this file.
00001 /* lu.c  -  PA = LU factorisation with pivoting
00002  * $Id: lu.c 1612 2006-06-01 19:04:35Z dangelo $
00003  * Copyright (C) 2000
00004  *    Antoine Lefebvre <antoine.lefebvre@polymtl.ca>
00005  *
00006  *
00007  * Licensed under the GPLv2
00008  */
00009    
00010 #include <stdlib.h>
00011 #include <stdio.h>
00012 #include <math.h>
00013 
00014 /*
00015   
00016    This algorithm will compute an LU factorisation on the augmented
00017    matrix (A) passed in arguments.
00018 
00019    This algorithm assumed the element on the diagonal of the lower
00020    triangular matrix set to 1.
00021   
00022    In order to save memory space, every coefficient of both matrix
00023    L and U are written in the original matrix overwriting initial
00024    values of A. 
00025 
00026 */
00027 
00028 int math_lu_solve(double *matrix, double *solution, int neq)
00029 {
00030   int i, j, k;
00031   
00032   int    idx;   /* index of the larger pivot */
00033   double big;       /* the larger pivot found */
00034   double tmp = 0.0;
00035   int itmp;
00036   
00037   int    *P;        /* keep memory of permutation (column permutation) */
00038   double *y;
00039     
00040   P = (int *) calloc (neq, sizeof(int));
00041   y = (double *) calloc (neq, sizeof(double));
00042 
00043   for (i = 0; i < neq; i++)
00044   {
00045     solution[i]  = 0; /* reset the solution vector */
00046     P[i] = i;         /* initialize permutation vector */
00047   }
00048     
00049   /* LU Factorisation */
00050 
00051   for (i = 0; i < neq - 1; i++) /* line */
00052   {
00053     for (j = i; j < neq; j++) /* column */
00054     {  
00055       tmp = 0.0;
00056       for (k = 0; k < i; k++)
00057         tmp += matrix[i + neq*P[k]] * matrix[k + neq*P[j]];
00058 
00059       matrix[i + neq*P[j]] = matrix[i + neq*P[j]] - tmp; 
00060     }
00061 
00062     /* find the larger pivot and interchange the columns */
00063     big = 0.0;
00064     idx = i;
00065     for (j = i; j < neq; j++)
00066     {
00067       if (big < fabs(matrix[i + neq*P[j]])) /* we found a larger pivot */
00068       {
00069         idx = j;
00070         big = fabs(matrix[i + neq*P[j]]);
00071       }
00072     }
00073     /* check if we have to interchange the lines */
00074     if (idx != i)
00075     {
00076       itmp   = P[i];
00077       P[i]   = P[idx];
00078       P[idx] = itmp;
00079     }
00080 
00081     if (matrix[i + neq*P[i]] == 0.0)
00082     {
00083       //printf("LU: matrix is singular, no unique solution.\n");
00084       free (y);
00085       free (P);
00086       return 0;
00087     }
00088     
00089     for (j = i+1; j < neq; j++)
00090     {
00091       tmp = 0.0;
00092       for (k = 0; k < i; k++)
00093         tmp += matrix[j + neq*P[k]] * matrix[k + neq*P[i]];
00094 
00095       matrix[j + neq*P[i]] = (matrix[j + neq*P[i]] - tmp)/matrix[i + neq*P[i]];
00096     }
00097   }
00098 
00099   i = neq - 1;  
00100   tmp = 0.0;
00101   for (k = 0; k < neq-1; k++)
00102     tmp += matrix[i + neq*P[k]] * matrix[k + neq*P[i]];
00103 
00104   matrix[i + neq*P[i]] = matrix[i + neq*P[i]] - tmp;
00105 
00106   
00107   /* End LU-Factorisation */
00108     
00109   /* substitution  for y    Ly = b*/
00110   for (i = 0; i < neq; i++)
00111   {
00112     tmp = 0.0;
00113     for (j = 0; j < i; j++)
00114       tmp += matrix[i + neq*P[j]] * y[j];
00115     
00116     y[i] = matrix[i + neq*neq] - tmp;
00117   }
00118   
00119   /* substitution for x   Ux = y*/
00120   for (i = neq - 1; i >=0; i--)
00121   {
00122     if (matrix[i + neq*P[i]] == 0.0)
00123     {
00124       //printf("LU: No unique solution exist.\n");
00125       free (y);
00126       free (P);
00127       return 0;
00128     }
00129     
00130     tmp = 0.0;
00131     for (j = i; j < neq; j++)
00132       tmp += matrix[i + neq*P[j]] * solution[P[j]];
00133     
00134     solution[P[i]] = (y[i] - tmp)/matrix[i + neq*P[i]];    
00135   }
00136      
00137   free (y);
00138   free (P);
00139   return 1;      
00140 }
00141 
00142 /* This function will divide each row of the matrix
00143  * by the highest element of this row.
00144  * This kind of scaling could improve precision while
00145  * solving some difficult matrix.
00146  */
00147 int NUM_matscale(double *matrix, int neq)
00148 {
00149   int i; /* line */
00150   int j; /* column */
00151 
00152   double val;
00153   double tmp;
00154   
00155   for (i = 0; i < neq; i++)
00156   {
00157     val = 0;
00158     /* find the highest value */
00159     for (j = 0; j < neq; j++)
00160     {
00161       tmp = fabs(matrix[i + neq*j]);
00162       val = (tmp > val) ? tmp : val;
00163     }
00164 
00165     /* divide element of the line by this value
00166      * including the right side
00167      * if the max value is defferent than zero
00168      */
00169     if (val != 0.0)
00170     {
00171       for (j = 0; j < neq+1; j++)
00172       matrix[i + neq*j] = matrix[i + neq*j]/val;
00173     }
00174   }
00175   return 0;
00176 
00177 }

Generated on Thu Jul 31 01:25:41 2014 for Hugintrunk by  doxygen 1.3.9.1