TABLE OF CONTENTS

ABINIT/lubksb [ Functions ]

[ Top ] [ Functions ]

NAME

  lubksb

FUNCTION

COPYRIGHT

  Copyright (C) 2014-2017 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

156 SUBROUTINE lubksb(a,n,np,indx,b)
157 
158  use m_profiling_abi
159 
160 !This section has been created automatically by the script Abilint (TD).
161 !Do not modify the following lines by hand.
162 #undef ABI_FUNC
163 #define ABI_FUNC 'lubksb'
164 !End of the abilint section
165 
166       IMPLICIT NONE
167 
168       INTEGER n,np,indx(n)
169       REAL*8 a(np,np),b(n)
170 
171 !     Solves the set of n linear equations A . X = B. Here a is input, not as
172 !     the matrix A but rather as its LU decomposition, determined by the
173 !     routine ludcmp. indx is input as the permutation vector returned by
174 !     ludcmp. b(1:n) is input as the right-hand side vector B, and returns
175 !     with the solution vector X. a, n, np, and indx are not modified by this
176 !     routine and can be left in place for successive calls with different
177 !     right-hand sides b. This routine takes into account the possibility that
178 !     b will begin with many zero elements, so it is efficient for use in
179 !     matrix inversion.
180 
181       INTEGER i,ii,j,ll
182       REAL*8 sum
183 !      write(std_out,*) 'ENTERING LUBKSB...'
184 !      write(std_out,201) ((a(i,j),j=1,n),i=1,n)
185 ! 201  FORMAT('A in lubksb ',/,3F16.8,/,3F16.8,/,3F16.8)
186 
187       ii=0
188       do i=1,n
189         ll=indx(i)
190         sum=b(ll)
191         b(ll)=b(i)
192         if (ii.ne.0)then
193           do j=ii,i-1
194             sum=sum-a(i,j)*b(j)
195           enddo 
196         else if (sum.ne.0.) then
197           ii=i 
198         endif
199         b(i)=sum
200       enddo 
201       do i=n,1,-1
202         sum=b(i)
203         do j=i+1,n
204           sum=sum-a(i,j)*b(j)
205         enddo
206         b(i)=sum/a(i,i) 
207       enddo
208 !      write(std_out,*) 'LEAVING LUBKSB...'
209       return 
210 END SUBROUTINE LUBKSB

ABINIT/ludcmp [ Functions ]

[ Top ] [ Functions ]

NAME

  ludcmp

FUNCTION

COPYRIGHT

  Copyright (C) 2014-2017 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

 28 #if defined HAVE_CONFIG_H
 29 #include "config.h"
 30 #endif
 31 
 32 #include "abi_common.h"
 33 
 34 
 35 SUBROUTINE ludcmp(a,n,np,indx,id,info)
 36 
 37  use defs_basis, only : std_out
 38 
 39 !This section has been created automatically by the script Abilint (TD).
 40 !Do not modify the following lines by hand.
 41 #undef ABI_FUNC
 42 #define ABI_FUNC 'ludcmp'
 43 !End of the abilint section
 44 
 45       implicit none
 46 
 47       INTEGER n,np,indx(n),NMAX,id,info
 48       REAL*8 a(np,np),TINY
 49       PARAMETER (NMAX=500,TINY=1.0e-20)
 50 
 51 !     Given a matrix a(1:n,1:n), with physical dimension np by np, this
 52 !     routine replaces it by the LU decomposition of a rowwise permutation of
 53 !     itself. a and n are input. a is output, arranged as in equation (2.3.14)
 54 !     above; indx(1:n) is an output vector that records the row permutation
 55 !     effected by the partial pivoting; id is output as +- 1 depending on
 56 !     whether the number of row interchanges was even or odd,
 57 !     respectively. This routine is used in combination with lubksb to solve
 58 !     linear equations or invert a matrix.
 59 
 60       INTEGER i,imax,j,k
 61       REAL*8 aamax,dum,sum,vv(NMAX) 
 62 
 63 !      write(std_out,*) 'ENTERING LUDCMP...'
 64 !      write(std_out,*) 'in ludcmp n=',n,' np=',np
 65 !      write(std_out,201) ((a(i,j),j=1,n),i=1,n)
 66 ! 201  FORMAT('A in ludcmp ',/,3F16.8,/,3F16.8,/,3F16.8)
 67       id=1
 68       info=0
 69       do i=1,n 
 70         aamax=0.
 71         do j=1,n
 72           if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
 73         enddo
 74         if (aamax.eq.0.) then
 75           write(std_out,*) 'LUDCMP: singular matrix !!!'
 76           do j=1,3
 77             write(std_out,*) (a(j,k),k=1,3)
 78           enddo
 79           info=1
 80           return
 81 !          stop 'singular matrix in ludcmp' 
 82         endif
 83         vv(i)=1./aamax 
 84       enddo 
 85       do j=1,n 
 86         do i=1,j-1 
 87           sum=a(i,j)
 88           do k=1,i-1
 89             sum=sum-a(i,k)*a(k,j)
 90           enddo 
 91           a(i,j)=sum
 92         enddo 
 93         aamax=0. 
 94         do i=j,n 
 95           sum=a(i,j)
 96           do k=1,j-1
 97             sum=sum-a(i,k)*a(k,j)
 98           enddo 
 99           a(i,j)=sum
100           dum=vv(i)*abs(sum) 
101           if (dum.ge.aamax) then 
102             imax=i
103             aamax=dum
104           endif
105         enddo 
106         if (j.ne.imax)then 
107           do  k=1,n 
108             dum=a(imax,k)
109             a(imax,k)=a(j,k)
110             a(j,k)=dum
111           enddo 
112           id=-id 
113           vv(imax)=vv(j) 
114         endif
115         indx(j)=imax
116         if(a(j,j).eq.0.)a(j,j)=TINY
117         if(j.ne.n)then 
118           dum=1./a(j,j)
119           do i=j+1,n
120             a(i,j)=a(i,j)*dum
121           enddo
122         endif
123       enddo 
124 !      write(std_out,*) 'LEAVING LUDCMP...'
125       return
126 END