1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2! Copyright 2010.  Los Alamos National Security, LLC. This material was    !
3! produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos !
4! National Laboratory (LANL), which is operated by Los Alamos National     !
5! Security, LLC for the U.S. Department of Energy. The U.S. Government has !
6! rights to use, reproduce, and distribute this software.  NEITHER THE     !
7! GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY,     !
8! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS         !
9! SOFTWARE.  If software is modified to produce derivative works, such     !
10! modified software should be clearly marked, so as not to confuse it      !
11! with the version available from LANL.                                    !
12!                                                                          !
13! Additionally, this program is free software; you can redistribute it     !
14! and/or modify it under the terms of the GNU General Public License as    !
15! published by the Free Software Foundation; version 2.0 of the License.   !
16! Accordingly, this program is distributed in the hope that it will be     !
17! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of   !
18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General !
19! Public License for more details.                                         !
20!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21
22SUBROUTINE SOLVEMATCG
23
24  USE CONSTANTS_MOD
25  USE FERMICOMMON
26  USE SETUPARRAY
27  USE SPINARRAY
28  USE MYPRECISION
29
30  IMPLICIT NONE
31
32  INTEGER :: I, J, ITER, BREAKLOOP
33  REAL(LATTEPREC) :: ERROR2
34  REAL(LATTEPREC) :: R0VEC, P0VEC, R1VEC, XALPHA, XBETA
35
36  IF (SPINON .EQ. 0) THEN
37
38#ifdef DOUBLEPREC
39     CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
40          BO, HDIM, BO, HDIM, 0.0D0, X2, HDIM)
41#elif defined(SINGLEPREC)
42     CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
43          BO, HDIM, BO, HDIM, 0.0, X2, HDIM)
44#endif
45
46     A = TWO*(X2 - BO)
47
48     DO I = 1, HDIM
49        A(I,I) = A(I,I) + ONE
50     ENDDO
51
52#ifdef DOUBLEPREC
53     CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
54          A, HDIM, BO, HDIM, 0.0D0, TMPMAT, HDIM)
55#elif defined(SINGLEPREC)
56     CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
57          A, HDIM, BO, HDIM, 0.0, TMPMAT, HDIM)
58#endif
59
60     R0 = TMPMAT - X2
61     P0 = MINUSONE*R0
62
63     ITER = 0
64
65     BREAKLOOP = 0
66
67     DO WHILE (BREAKLOOP .EQ. 0)
68
69        ITER = ITER + 1
70
71        IF (ITER .EQ. 50) THEN
72           CALL PANIC
73           CALL ERRORS("solvematcg_sparse",'("SOLVEMATCG NOT CONVERGING")')
74        ENDIF
75
76
77#ifdef DOUBLEPREC
78        CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
79             A, HDIM, P0, HDIM, 0.0D0, TMPMAT, HDIM)
80#elif defined(SINGLEPREC)
81        CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
82             A, HDIM, P0, HDIM, 0.0, TMPMAT, HDIM)
83#endif
84
85        ERROR2 = ZERO
86
87        !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(GUIDED) &
88        !$OMP SHARED (TMPMAT, R0, P0, BO,HDIM) &
89        !$OMP PRIVATE(I, J, R0VEC, P0VEC, R1VEC, XALPHA, XBETA) &
90        !$OMP REDUCTION(+ : ERROR2)
91
92        DO I = 1, HDIM
93
94           R0VEC = ZERO
95           P0VEC = ZERO
96           R1VEC = ZERO
97
98           DO J = 1, HDIM
99
100              P0VEC = P0VEC + P0(J,I)*TMPMAT(J,I)
101              R0VEC = R0VEC + R0(J,I)*R0(J,I)
102
103           ENDDO
104
105           XALPHA = R0VEC/P0VEC
106
107           DO J = 1, HDIM
108
109              ! New density matrix
110
111              BO(J,I) = BO(J,I) + P0(J,I)*XALPHA
112
113              ! Calculating R1
114
115              R0(J,I) = R0(J,I) + TMPMAT(J,I)*XALPHA
116
117              R1VEC = R1VEC + R0(J,I)*R0(J,I)
118
119           ENDDO
120
121           ERROR2 = ERROR2 + R1VEC
122
123           XBETA = R1VEC/R0VEC
124
125           DO J = 1, HDIM
126
127              P0(J,I) = P0(J,I)*XBETA - R0(J,I)
128
129           ENDDO
130
131        ENDDO
132
133        !$OMP END PARALLEL DO
134
135        IF (ERROR2 .LT. CGTOL2) THEN
136           BREAKLOOP = 1
137        ENDIF
138
139        !        PRINT*, ITER, ERROR2
140
141     ENDDO
142
143  ELSE
144
145     ! This is the spin-polarized version
146
147#ifdef DOUBLEPREC
148     CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
149          RHOUP, HDIM, RHOUP, HDIM, 0.0D0, X2, HDIM)
150#elif defined(SINGLEPREC)
151     CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
152          RHOUP, HDIM, RHOUP, HDIM, 0.0, X2, HDIM)
153#endif
154
155     A = TWO*(X2 - RHOUP)
156
157     DO I = 1, HDIM
158        A(I,I) = A(I,I) + ONE
159     ENDDO
160
161#ifdef DOUBLEPREC
162     CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
163          A, HDIM, RHOUP, HDIM, 0.0D0, TMPMAT, HDIM)
164#elif defined(SINGLEPREC)
165     CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
166          A, HDIM, RHOUP, HDIM, 0.0, TMPMAT, HDIM)
167#endif
168
169     R0 = TMPMAT - X2
170     P0 = MINUSONE*R0
171
172     ITER = 0
173
174     BREAKLOOP = 0
175
176     DO WHILE (BREAKLOOP .EQ. 0)
177
178        ITER = ITER + 1
179
180        IF (ITER .EQ. 50) THEN
181           CALL PANIC
182           CALL ERRORS("solvematcg_sparse",'("SOLVEMATCG NOT CONVERGING: SPIN UP")')
183        ENDIF
184
185
186        !     PRINT*, ITER, ERROR2
187
188#ifdef DOUBLEPREC
189        CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
190             A, HDIM, P0, HDIM, 0.0D0, TMPMAT, HDIM)
191#elif defined(SINGLEPREC)
192        CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
193             A, HDIM, P0, HDIM, 0.0, TMPMAT, HDIM)
194#endif
195
196        ERROR2 = ZERO
197
198        !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(GUIDED) &
199        !$OMP SHARED (TMPMAT, R0, P0, RHOUP, HDIM) &
200        !$OMP PRIVATE(I, J, R0VEC, P0VEC, R1VEC, XALPHA, XBETA) &
201        !$OMP REDUCTION(+ : ERROR2)
202
203        DO I = 1, HDIM
204
205           R0VEC = ZERO
206           P0VEC = ZERO
207           R1VEC = ZERO
208
209           DO J = 1, HDIM
210
211              P0VEC = P0VEC + P0(J,I)*TMPMAT(J,I)
212              R0VEC = R0VEC + R0(J,I)*R0(J,I)
213
214           ENDDO
215
216           XALPHA = R0VEC/P0VEC
217
218           DO J = 1, HDIM
219
220              ! New density matrix
221
222              RHOUP(J,I) = RHOUP(J,I) + P0(J,I)*XALPHA
223
224              ! Calculating R1
225
226              R0(J,I) = R0(J,I) + TMPMAT(J,I)*XALPHA
227
228              R1VEC = R1VEC + R0(J,I)*R0(J,I)
229
230           ENDDO
231
232           ERROR2 = ERROR2 + R1VEC
233
234           XBETA = R1VEC/R0VEC
235
236           DO J = 1, HDIM
237
238              P0(J,I) = P0(J,I)*XBETA - R0(J,I)
239
240           ENDDO
241
242        ENDDO
243
244        !$OMP END PARALLEL DO
245
246        !        PRINT*, "UP ", ITER, ERROR2
247
248        IF (ITER .GT. 3 .AND. ERROR2 .LT. CGTOL2) THEN
249           BREAKLOOP = 1
250        ENDIF
251
252     ENDDO
253
254#ifdef DOUBLEPREC
255     CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
256          RHODOWN, HDIM, RHODOWN, HDIM, 0.0D0, X2, HDIM)
257#elif defined(SINGLEPREC)
258     CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
259          RHODOWN, HDIM, RHODOWN, HDIM, 0.0, X2, HDIM)
260#endif
261
262     A = TWO*(X2 - RHODOWN)
263
264     DO I = 1, HDIM
265        A(I,I) = A(I,I) + ONE
266     ENDDO
267
268#ifdef DOUBLEPREC
269     CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
270          A, HDIM, RHODOWN, HDIM, 0.0D0, TMPMAT, HDIM)
271#elif defined(SINGLEPREC)
272     CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
273          A, HDIM, RHODOWN, HDIM, 0.0, TMPMAT, HDIM)
274#endif
275
276     R0 = TMPMAT - X2
277     P0 = MINUSONE*R0
278
279     ITER = 0
280
281     BREAKLOOP = 0
282
283     DO WHILE  (BREAKLOOP .EQ. 0)
284
285        ITER = ITER + 1
286
287        IF (ITER .EQ. 50) THEN
288           CALL PANIC
289           CALL ERRORS("solvematcg_sparse",'("SOLVEMATCG NOT CONVERGING: SPIN DOWN")')
290        ENDIF
291
292        !     PRINT*, ITER, ERROR2
293
294#ifdef DOUBLEPREC
295        CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0D0, &
296             A, HDIM, P0, HDIM, 0.0D0, TMPMAT, HDIM)
297#elif defined(SINGLEPREC)
298        CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, 1.0, &
299             A, HDIM, P0, HDIM, 0.0, TMPMAT, HDIM)
300#endif
301
302        ERROR2 = ZERO
303
304        !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(GUIDED) &
305        !$OMP SHARED (TMPMAT, R0, P0, RHODOWN, HDIM) &
306        !$OMP PRIVATE(I, J, R0VEC, P0VEC, R1VEC, XALPHA, XBETA) &
307        !$OMP REDUCTION(+ : ERROR2)
308
309        DO I = 1, HDIM
310
311           R0VEC = ZERO
312           P0VEC = ZERO
313           R1VEC = ZERO
314
315           DO J = 1, HDIM
316
317              P0VEC = P0VEC + P0(J,I)*TMPMAT(J,I)
318              R0VEC = R0VEC + R0(J,I)*R0(J,I)
319
320           ENDDO
321
322           XALPHA = R0VEC/P0VEC
323
324           DO J = 1, HDIM
325
326              ! New density matrix
327
328              RHODOWN(J,I) = RHODOWN(J,I) + P0(J,I)*XALPHA
329
330              ! Calculating R1
331
332              R0(J,I) = R0(J,I) + TMPMAT(J,I)*XALPHA
333
334              R1VEC = R1VEC + R0(J,I)*R0(J,I)
335
336           ENDDO
337
338           ERROR2 = ERROR2 + R1VEC
339
340           XBETA = R1VEC/R0VEC
341
342           DO J = 1, HDIM
343
344              P0(J,I) = P0(J,I)*XBETA - R0(J,I)
345
346           ENDDO
347
348        ENDDO
349
350        !$OMP END PARALLEL DO
351
352        !        PRINT*, "DOWN ", ITER, ERROR2
353
354        IF (ITER .GT. 3 .AND. ERROR2 .LT. CGTOL2) THEN
355           BREAKLOOP = 1
356        ENDIF
357
358     ENDDO
359
360  ENDIF
361
362  RETURN
363
364END SUBROUTINE SOLVEMATCG
365