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