Chimie Théorique » Opt'n Path

      SUBROUTINE gaussj(a,n,np,b,m,mp)

! This routine solves the system AX=B by inverting A.

! This is based on Gauss-Jordan pivoting.

! Input:

! A is NxN REAL(KREAL) matrix, order N

! B is NxM REAL(KREAL)  matrix (each ocolumn j of this B matrix corresponds to

!   a AX=B_j problem in fact).

! np: actual size of A(np,np)

! mp: actual size of B(np,mp)

!

! Output:

! A is destroyed and contains A inverse

! B is destroyed and contains the solution X vectors.

!

! Adapted from  "Numerical Recipes in F77"

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs 

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" 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 Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

     IMPLICIT NONE

     INTEGER, PARAMETER :: KINT=KIND(1)

     INTEGER, PARAMETER :: KREAL=KIND(1.0D0)

      INTEGER(KINT) :: m,mp,n,np,NMAX

      REAL(KREAL) :: a(np,np),b(np,mp)

      PARAMETER (NMAX=50)

      INTEGER(KINT) :: i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)

      REAL(KREAL) :: big,dum,pivinv

      do 11 j=1,n

        ipiv(j)=0

11    continue

      do 22 i=1,n

        big=0.d0

        do 13 j=1,n

          if(ipiv(j).ne.1)then

            do 12 k=1,n

              if (ipiv(k).eq.0) then

                if (abs(a(j,k)).ge.big)then

                  big=abs(a(j,k))

                  irow=j

                  icol=k

                endif

              else if (ipiv(k).gt.1) then

                Write(*,*) 'singular matrix in gaussj'

                STOP

              endif

12          continue

          endif

13      continue

        ipiv(icol)=ipiv(icol)+1

        if (irow.ne.icol) then

          do 14 l=1,n

            dum=a(irow,l)

            a(irow,l)=a(icol,l)

            a(icol,l)=dum

14        continue

          do 15 l=1,m

            dum=b(irow,l)

            b(irow,l)=b(icol,l)

            b(icol,l)=dum

15        continue

        endif

        indxr(i)=irow

        indxc(i)=icol

        if (a(icol,icol).eq.0.d0) THEN

           WRITE(*,*)  'singular matrix in gaussj'

           STOP

        END IF

        pivinv=1.d0/a(icol,icol)

        a(icol,icol)=1.d0

        do 16 l=1,n

          a(icol,l)=a(icol,l)*pivinv

16      continue

        do 17 l=1,m

          b(icol,l)=b(icol,l)*pivinv

17      continue

        do 21 ll=1,n

          if(ll.ne.icol)then

            dum=a(ll,icol)

            a(ll,icol)=0.d0

            do 18 l=1,n

              a(ll,l)=a(ll,l)-a(icol,l)*dum

18          continue

            do 19 l=1,m

              b(ll,l)=b(ll,l)-b(icol,l)*dum

19          continue

          endif

21      continue

22    continue

      do 24 l=n,1,-1

        if(indxr(l).ne.indxc(l))then

          do 23 k=1,n

            dum=a(k,indxr(l))

            a(k,indxr(l))=a(k,indxc(l))

            a(k,indxc(l))=dum

23        continue

        endif

24    continue

      return

      END