1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8
9      SUBROUTINE INVER(A,B,N,NDIM,INFO)
10
11      IMPLICIT NONE
12      INTEGER N,NDIM,INFO
13      DOUBLE PRECISION A(NDIM,NDIM),B(NDIM,NDIM),X
14
15C Routine to obtain the inverse of a general, nonsymmetric
16C square matrix, by calling the appropriate LAPACK routines
17C The matrix A is of order N, but is defined with a
18C size NDIM >= N.
19C If the LAPACK routines fail, try the good old Numerical Recipes
20C algorithm
21C
22C P. Ordejon, June 2003
23C
24C **** INPUT ****
25C A(NDIM,NDIM)   real*8      Matrix to be inverted
26C N              integer     Size of the matrix
27C NDIM           integer     Defined size of matrix
28C **** OUTPUT ****
29C B(NDIM,NDIM)   real*8      Inverse of A
30C INFO           integer     If inversion was unsucessful,
31C                            INFO <> 0 is returned
32C                            Is successfull, INFO = 0 returned
33C ***************
34
35
36      DOUBLE PRECISION WORK(N),C,ERROR,DELTA,TOL
37      INTEGER IPIV(N)
38      INTEGER I,J,K
39
40      TOL=1.0D-4
41      INFO = 0
42
43      DO I=1,N
44      DO J=1,N
45        B(I,J)=A(I,J)
46      ENDDO
47      ENDDO
48
49      CALL DGETRF(N,N,B,NDIM,IPIV,INFO)
50
51      IF (INFO .NE. 0) THEN
52C        'inver:  ERROR: DGETRF exited with error message', INFO
53        GOTO 100
54      ENDIF
55
56      CALL DGETRI(N,B,NDIM,IPIV,WORK,N,INFO)
57
58      IF (INFO .NE. 0) THEN
59C     .    'inver:  ERROR: DGETRI exited with error message', INFO
60        GOTO 100
61      ENDIF
62
63C CHECK THAT THE INVERSE WAS CORRECTLY CALCULATED
64
65
66      ERROR=0.0D0
67      DO I=1,N
68      DO J=1,N
69        C=0.0D0
70        DO K=1,N
71          C=C+A(I,K)*B(K,J)
72        ENDDO
73        DELTA=0.0D0
74        IF (I.EQ.J)  DELTA=1.0D0
75        ERROR=ERROR+DABS(C-DELTA)
76        C=0.0D0
77        DO K=1,N
78          C=C+B(I,K)*A(K,J)
79        ENDDO
80        DELTA=0.0D0
81        IF (I.EQ.J)  DELTA=1.0D0
82        ERROR=ERROR+DABS(C-DELTA)
83      ENDDO
84      ENDDO
85
86      IF (ERROR/N .GT. TOL) THEN
87        GOTO 100
88      ENDIF
89
90      INFO = 0
91      RETURN
92
93100   CONTINUE
94
95C Try simple, old algorithm:
96
97      DO I=1,N
98        DO J=1,N
99          B(I,J)=A(I,J)
100        ENDDO
101      ENDDO
102      DO I=1,N
103        IF (B(I,I) .EQ. 0.0D0) THEN
104          INFO = 1
105          RETURN
106        ENDIF
107        X=B(I,I)
108        B(I,I)=1.0d0
109        DO J=1,N
110          B(J,I)=B(J,I)/X
111        ENDDO
112        DO K=1,N
113          IF ((K-I).NE.0) THEN
114            X=B(I,K)
115            B(I,K)=0.0d0
116            DO J=1,N
117              B(J,K)=B(J,K)-B(J,I)*X
118            ENDDO
119          ENDIF
120        ENDDO
121      ENDDO
122
123C CHECK THAT THE INVERSE WAS CORRECTLY CALCULATED
124
125      ERROR=0.0D0
126      DO I=1,N
127      DO J=1,N
128        C=0.0D0
129        DO K=1,N
130          C=C+A(I,K)*B(K,J)
131        ENDDO
132        DELTA=0.0D0
133        IF (I.EQ.J)  DELTA=1.0D0
134        ERROR=ERROR+DABS(C-DELTA)
135        C=0.0D0
136        DO K=1,N
137          C=C+B(I,K)*A(K,J)
138        ENDDO
139        DELTA=0.0D0
140        IF (I.EQ.J)  DELTA=1.0D0
141        ERROR=ERROR+DABS(C-DELTA)
142      ENDDO
143      ENDDO
144
145      IF (ERROR/N .GT. TOL) THEN
146C        WRITE(6,'(a,f10.5)')
147C     .    'inver:  INVER unsuccessful, ERROR = ', ERROR
148        INFO = 1
149        RETURN
150      ENDIF
151
152
153
154      INFO = 0
155      RETURN
156      END
157