[cig-commits] r6178 - short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d

[willic3 at geodynamics.org]

Author: willic3
Date: 2007-03-05 18:32:51 -0800 (Mon, 05 Mar 2007)
New Revision: 6178

Added:
   short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f
   short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f
Modified:
   short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am
Log:
Temporarily added Numerical Recipes-derived routines for LU
decomposition and backsubstitution.  This should be replaced with
LAPACK calls to routines that take advantage of matrix symmetry and
packed storage.



Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am	2007-03-05 18:55:00 UTC (rev 6177)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am	2007-03-06 02:32:51 UTC (rev 6178)
@@ -98,6 +98,8 @@
 	local.f \
 	localf.f \
 	localx.f \
+	lubksb.f \
+	ludcmp.f \
 	makemsr.F \
 	mat_1.f \
 	mat_2.f \

Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f	2007-03-05 18:55:00 UTC (rev 6177)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f	2007-03-06 02:32:51 UTC (rev 6178)
@@ -0,0 +1,41 @@
+      SUBROUTINE lubksb(a,n,np,indx,b)
+c
+c...  routine to perform backsubstitution on an LU-decomposed matrix.
+c     Adapted from Numerical Recipes.
+c
+      include "implicit.inc"
+c
+c...  subroutine arguments
+c
+      integer n,np
+      INTEGER indx(n)
+      double precision a(np,np),b(n)
+c
+c...  local variables
+c
+      INTEGER i,ii,j,ll
+      double precision sum
+c
+      ii=0
+      do 12 i=1,n
+        ll=indx(i)
+        sum=b(ll)
+        b(ll)=b(i)
+        if (ii.ne.0)then
+          do 11 j=ii,i-1
+            sum=sum-a(i,j)*b(j)
+11        continue
+        else if (sum.ne.0.0d0) then
+          ii=i
+        endif
+        b(i)=sum
+12    continue
+      do 14 i=n,1,-1
+        sum=b(i)
+        do 13 j=i+1,n
+          sum=sum-a(i,j)*b(j)
+13      continue
+        b(i)=sum/a(i,i)
+14    continue
+      return
+      END

Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f	2007-03-05 18:55:00 UTC (rev 6177)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f	2007-03-06 02:32:51 UTC (rev 6178)
@@ -0,0 +1,73 @@
+      SUBROUTINE ludcmp(a,n,np,indx,d)
+c
+c...  subroutine to perform an LU decomposition on a matrix.
+c     Adapted from Numerical Recipes.
+c
+      include "implicit.inc"
+c
+c...  parameter definitions
+c
+      integer NMAX
+      double precision TINY
+      PARAMETER (NMAX=500,TINY=1.0d-20)
+c
+c...  subroutine arguments
+c
+      integer n,np
+      INTEGER indx(n)
+      double precision d,a(np,np)
+c
+c...  local variables
+c
+      INTEGER i,imax,j,k
+      double precision aamax,dum,sum,vv(NMAX)
+      d=1.0d0
+      do 12 i=1,n
+        aamax=0.0d0
+        do 11 j=1,n
+          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
+11      continue
+        if (aamax.eq.0.0d0) pause 'singular matrix in ludcmp'
+        vv(i)=1.0d0/aamax
+12    continue
+      do 19 j=1,n
+        do 14 i=1,j-1
+          sum=a(i,j)
+          do 13 k=1,i-1
+            sum=sum-a(i,k)*a(k,j)
+13        continue
+          a(i,j)=sum
+14      continue
+        aamax=0.0d0
+        do 16 i=j,n
+          sum=a(i,j)
+          do 15 k=1,j-1
+            sum=sum-a(i,k)*a(k,j)
+15        continue
+          a(i,j)=sum
+          dum=vv(i)*abs(sum)
+          if (dum.ge.aamax) then
+            imax=i
+            aamax=dum
+          endif
+16      continue
+        if (j.ne.imax)then
+          do 17 k=1,n
+            dum=a(imax,k)
+            a(imax,k)=a(j,k)
+            a(j,k)=dum
+17        continue
+          d=-d
+          vv(imax)=vv(j)
+        endif
+        indx(j)=imax
+        if(a(j,j).eq.0.0d0)a(j,j)=TINY
+        if(j.ne.n)then
+          dum=1.0d0/a(j,j)
+          do 18 i=j+1,n
+            a(i,j)=a(i,j)*dum
+18        continue
+        endif
+19    continue
+      return
+      END