ludcmp.h

Go to the documentation of this file.

/*
   This file belongs to Aeneas. Aeneas is a GNU package released under GPL 3.
   This code is a simulator for Submicron 3D Semiconductor Devices. 
   It implements the Monte Carlo transport in 3D tetrahedra meshes
   for the simulation of the semiclassical Boltzmann equation for both electrons.
   It also includes all the relevant quantum effects for nanodevices.

   Copyright (C) 2007 Jean Michel Sellier <sellier@dmi.unict.it>
 
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 3, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

// Created on : 08 june 2007, Siracusa, Jean Michel Sellier
// Last modified : 16 september 2007, Siracusa, Jean Michel Sellier

// given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition
// of a rowwise permutation of itself. a and n are input. a is output, arranged
// as a lu decomposition; indx[1..n] is an output vector that records the row
// permutation effected by the partial pivoting; d is the output as +-1 depending
// on whether the number of row interchanges was even or odd, respectively.
// This routine is used in combination with lubksb to solve linear equations
// or invert a matrix.

{
 int i, imax, j, k;
 double big, dum, sum, temp;
 double *vv;
 
// imax=0; // dummy value, just to avoid warnings during the compilation...

 vv = malloc((n+1)*sizeof(double)); // vv[1..n]
 if(vv==NULL){
   printf("ludcmp error : not enough memory!\n");
   exit(0);
 }
/*     {int ii,jj;
     for(ii=1;ii<=4;ii++){
        for(jj=1;jj<=4;jj++) printf("*[%d][%d] = %g   ",ii,jj,a[ii][jj]);
        printf("\n");
     }}*/

 d=1.0;
 for(i=1;i<=n;i++){
    big=0.0;
    for(j=1;j<=n;j++)
      if((temp=fabs(a[i][j]))>big) big=temp;
    if(big==0.0){
     int ii,jj;
     for(ii=1;ii<=4;ii++){
        for(jj=1;jj<=4;jj++) printf("a[%d][%d] = %g   ",ii,jj,a[ii][jj]);
        printf("\n");
     }
     printf("ludcmp error : singular matrix!\n");
     exit(0);
    }
// non zero largest element
    vv[i]=1.0/big; // save the scaling    
 }
 for(j=1;j<=n;j++){
    for(i=1;i<j;i++){
      sum=a[i][j];
      for(k=1;k<i;k++) sum-=a[i][k]*a[k][j];
      a[i][j]=sum;
    }
    big=0.0; // initialise for the search for largest pivot element
    for(i=j;i<=n;i++){
      sum=a[i][j];
      for(k=1;k<j;k++) sum-=a[i][k]*a[k][j];
      a[i][j]=sum;
      if((dum=vv[i]*fabs(sum))>=big){
// Is the figure of merit for the pivot better than the best so far?
        big=dum;
        imax=i;
      }
    }
    if(j!=imax){ // do we need to interchage the row?
      for(k=1;k<=n;k++){
        dum=a[imax][k];
        a[imax][k]=a[j][k];
        a[j][k]=dum;
      }
      d=-d; // ...and change the parity of d.
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if(a[j][j]==0.0) a[j][j]=TINY;
// if the pivot element is zero the matrix is singular (at least to the precision
// of the algorithm). For some applications on singular matrices, it is desirable
// to substitute TINY for zero.
    if(j!=n){
      dum=1.0/(a[j][j]);
      for(i=j+1;i<=n;i++) a[i][j]*=dum;
    }
 }
 free(vv);
}

---

© sourcejam.com 2005-2008