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