1      SUBROUTINE DIIS_SIMPLE(LUEVEC,NVEC,NDIM,C)
2*
3* NVEC error-vectors are given on LUEVEC
4* Use DIIS to find optimal combination of these vectors
5*
6*. Simple and generic DIIS code - more simple than efficient
7*. All subspace matrices constructed without reuse
8*
9*. Jeppe Olsen, August 2005
10*
11      INCLUDE 'wrkspc.inc'
12*. Output : Optimial combination of directions
13       REAL*8 INPROD
14       DIMENSION C(NVEC)
15*. Non-redundant formalism : X = X(NVEC) + sum(i=1,NVEC-1)C(I)(X(I)-X(N))
16*  corresponding error       E = E(NVEC) + sum(i=1,NVEC-1)C(I)(E(I)-E(N))
17*
18* Minimizing norm of gradient corresponds to solving linear set of equations
19*
20* A C = RHS
21* A_ij = <E(I)-E(N)!E(J)-E(N)>, RHS_i = -<E(I)-E(N)!E(N)>
22*
23*
24      NTEST = 100
25*
26      IF(NTEST.GE.10) THEN
27        WRITE(6,*) ' Entering DIIS '
28        WRITE(6,*) ' Number of vectors ', NVEC
29      END IF
30*
31      IDUM = 0
32      CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'DIISS')
33*. A few scratch vectors
34      LEN = MAX(NDIM,NVEC-1)
35      CALL MEMMAN(KLVEC1,LEN,'ADDL ',2,'VC1_DI')
36      CALL MEMMAN(KLVEC2,NDIM,'ADDL ',2,'VC2_DI')
37      CALL MEMMAN(KLVEC3,NDIM,'ADDL ',2,'VC3_DI')
38*. Space for A and RHS and for inverting
39      CALL MEMMAN(KLA,(NVEC-1)**2,'ADDL ',2,'AMAT  ')
40      CALL MEMMAN(KLRHS,(NVEC-1), 'ADDL ',2,'RHS   ')
41      CALL MEMMAN(KLSCRM1,(NVEC-1)**2,'ADDL ',2,'SCRM1 ')
42      CALL MEMMAN(KLSCRM2,(NVEC-1)**2,'ADDL ',2,'SCRM2 ')
43
44*. Read E(NVEC) in VEC1
45      CALL REWINO(LUEVEC)
46      DO IVEC = 1, NVEC
47C?      WRITE(6,*) ' Initial loop, IVEC = ', IVEC
48C  VEC_FROM_DISC(VEC,LENGTH,IREW,LBLK,LU)
49        CALL VEC_FROM_DISC(WORK(KLVEC1),NDIM,0,-1,LUEVEC)
50      END DO
51*
52      ONE = 1.0D0
53      ONEM = -1.0D0
54      DO IVEC = 1, NVEC-1
55C?     WRITE(6,*) ' IVEC = ', IVEC
56*. Read vector E(IVEC) in VEC2
57       CALL REWINO(LUEVEC)
58       DO KVEC = 1, IVEC
59         CALL VEC_FROM_DISC(WORK(KLVEC2),NDIM,0,-1,LUEVEC)
60       END DO
61*. E(IVEC)-E(NVEC) in VEC2
62       CALL VECSUM(WORK(KLVEC2),WORK(KLVEC2),WORK(KLVEC1),
63     &             ONE,ONEM,NDIM)
64       WORK(KLRHS-1+IVEC) = INPROD(WORK(KLVEC2),WORK(KLVEC1),NDIM)
65       EEII = INPROD(WORK(KLVEC2),WORK(KLVEC2),NDIM)
66       WORK(KLA-1+(IVEC-1)*(NVEC-1)+IVEC) = EEII
67       DO JVEC = IVEC+1, NVEC-1
68C?     WRITE(6,*) ' JVEC = ', JVEC
69*. Read vector E(JVEC) on VEC3
70         CALL VEC_FROM_DISC(WORK(KLVEC3),NDIM,0,-1,LUEVEC)
71*. E(JVEC) - E(NVEC) in VEC3
72         CALL VECSUM(WORK(KLVEC3),WORK(KLVEC3),WORK(KLVEC1),
73     &               ONE,ONEM,NDIM)
74         EEIJ = INPROD(WORK(KLVEC2),WORK(KLVEC3),NDIM)
75         WORK(KLA-1+(IVEC-1)*(NVEC-1)+JVEC) = EEIJ
76         WORK(KLA-1+(JVEC-1)*(NVEC-1)+IVEC) = EEIJ
77       END DO
78      END DO
79*. Well, the RHS should contain a - so
80      CALL SCALVE(WORK(KLRHS),ONEM,NVEC-1)
81*
82      IF(NTEST.GE.100) THEN
83        WRITE(6,*) 'The <E(I)-E(N)!E(J)-E(N)> matrix in DIIS '
84        CALL WRTMAT(WORK(KLA),NVEC-1,NVEC-1,NVEC-1,NVEC-1)
85        WRITE(6,*) 'The <E(I)-E(N)!E(N)> vector in DIIS '
86        CALL WRTMAT(WORK(KLRHS),1,NVEC-1,1,NVEC-1)
87      END IF
88*. Solve linear equations
89*. Invert EE
90      IF(NVEC.GE.1) THEN
91      CALL INVERT_BY_DIAG2(WORK(KLA),WORK(KLSCRM1),WORK(KLSCRM2),
92     &                    WORK(KLVEC1),NVEC-1)
93      END IF
94      CALL MATVCB(WORK(KLA),WORK(KLRHS),C,NVEC-1,NVEC-1,0)
95*. IN C, we have coefficients for C(I),I.ne.NVEC, obtain C(NVEC)
96      C(NVEC) = 1.0D0-ELSUM(C,NVEC-1)
97*. Clean up : leave the file at end of vector NVEC
98      IF(NVEC.GE.2) THEN
99         CALL VEC_FROM_DISC(WORK(KLVEC2),NDIM,0,-1,LUEVEC)
100      END IF
101*
102      IF(NTEST.GE.10) THEN
103        WRITE(6,*) ' The C-coefficients obtained from DIIS '
104        CALL WRTMAT(C,1,NVEC,1,NVEC)
105      END IF
106*
107C     WRITE(6,*) ' Leaving DIIS '
108      CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'DIISS')
109      RETURN
110      END
111
112
113
114c $Id$
115