Chimie Théorique » Opt'n Path

      SUBROUTINE MINV(A,N,D)

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

! This routine was adapted from the public domain mopac6 diis.f

!  source file (c) 2009, Stewart Computational Chemistry.

!  <http://www.openmopac.net/Downloads/Downloads.html>

!

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

!  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), INTENT(IN) :: N

      REAL(KREAL), INTENT(INOUT) :: A(N*N)

      REAL(KREAL), INTENT(OUT) :: D

!**********************************************************************

!*

!*     INVERT A MATRIX USING GAUSS-JORDAN METHOD.  PART OF DIIS

!*     A - INPUT MATRIX (MUST BE A GENERAL MATRIX), DESTROYED IN

!*        COMPUTATION AND REPLACED BY RESULTANT INVERSE.

!*     N - ORDER OF MATRIX A

!*     D - RESULTANT DETERMINANT

!*

!**********************************************************************

      INTEGER(KINT), ALLOCATABLE :: L(:), M(:)

      INTEGER(KINT) :: I,J,K,NK,KK,KI,KJ,IZ,IJ,JP,JK,IK,JR,JQ,JI

      REAL(KREAL) :: BigA,Hold

      ALLOCATE(L(N),M(N))

!

!     SEARCH FOR LARGEST ELEMENT

!

      D=1.0D0

      NK=-N

      DO 180 K=1,N

         NK=NK+N

         L(K)=K

         M(K)=K

         KK=NK+K

         BIGA=A(KK)

         DO 20 J=K,N

            IZ=N*(J-1)

            DO 20 I=K,N

               IJ=IZ+I

               IF (ABS(BIGA).LT.ABS(A(IJ)))THEN

                  BIGA=A(IJ)

                  L(K)=I

                  M(K)=J

               ENDIF

20          CONTINUE

!

!     INTERCHANGE ROWS

!

         J=L(K)

         IF (J-K) 50,50,30

   30    KI=K-N

         DO 40 I=1,N

            KI=KI+N

            HOLD=-A(KI)

            JI=KI-K+J

            A(KI)=A(JI)

   40    A(JI)=HOLD

!

!     INTERCHANGE COLUMNS

!

   50    I=M(K)

         IF (I-K) 80,80,60

   60    JP=N*(I-1)

         DO 70 J=1,N

            JK=NK+J

            JI=JP+J

            HOLD=-A(JK)

            A(JK)=A(JI)

   70    A(JI)=HOLD

!C

!C     DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS

!C     CONTAINED IN BIGA)

!C

   80    IF (BIGA) 100,90,100

   90    D=0.0

         RETURN

  100    DO 120 I=1,N

            IF (I-K) 110,120,110

  110       IK=NK+I

            A(IK)=A(IK)/(-BIGA)

  120    CONTINUE

!C  REDUCE MATRIX

         DO 150 I=1,N

            IK=NK+I

            HOLD=A(IK)

            IJ=I-N

            DO 150 J=1,N

               IJ=IJ+N

               IF (I-K) 130,150,130

  130          IF (J-K) 140,150,140

  140          KJ=IJ-I+K

               A(IJ)=HOLD*A(KJ)+A(IJ)

  150    CONTINUE

!C

!C     DIVIDE ROW BY PIVOT

!C

         KJ=K-N

         DO 170 J=1,N

            KJ=KJ+N

            IF (J-K) 160,170,160

  160       A(KJ)=A(KJ)/BIGA

  170    CONTINUE

!C

!C     PRODUCT OF PIVOTS

!C

         D=MAX(-1.D25,MIN(1.D25,D))

         D=D*BIGA

!C

!C     REPLACE PIVOT BY RECIPROCAL

!C

         A(KK)=1.0/BIGA

  180 CONTINUE

!C

!C     FINAL ROW AND COLUMN INTERCHANGE

!C

      K=N

  190 K=(K-1)

      IF (K) 260,260,200

  200 I=L(K)

      IF (I-K) 230,230,210

  210 JQ=N*(K-1)

      JR=N*(I-1)

      DO 220 J=1,N

         JK=JQ+J

         HOLD=A(JK)

         JI=JR+J

         A(JK)=-A(JI)

  220 A(JI)=HOLD

  230 J=M(K)

      IF (J-K) 190,190,240

  240 KI=K-N

      DO 250 I=1,N

         KI=KI+N

         HOLD=A(KI)

         JI=KI-K+J

         A(KI)=-A(JI)

  250 A(JI) =HOLD

      GO TO 190

  260 RETURN

      END