1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck ciham */ 20 SUBROUTINE CIHAM(IWAY,MXPRDT,ICSYM,IDET,CIDIA,NCONF, 21 * ECORE,FCAC,H2AC,XNDXCI,WRK,LWRK,IPRINT) 22C 23C 29-Jan-1989 hjaaj 24C 25C Explicit CI Hamiltonian. 26C 27C Notes: 28C - introduce symtest using IHSYM=0 29C 30#include "implicit.h" 31#include "priunit.h" 32 DIMENSION IDET(*),CIDIA(*),FCAC(*),H2AC(*),XNDXCI(*),WRK(LWRK) 33C 34C Used from common blocks: 35C INFINP : MCTYPE 36C INFORB : MULD2H(8,8),NASHT,N2ASHX 37C CIINFO : ICOMBI, IPSIGN 38C DETBAS : KLTSOB 39C 40#include "maxorb.h" 41#include "infinp.h" 42#include "inforb.h" 43#include "ciinfo.h" 44#include "detbas.h" 45C 46C Set some codes: 47C 48 PSIGN = IPSIGN 49 NPRDET = 0 50 IHSYM = 1 51C ... Hamiltonian matrix is always symmetric in this version. 52 IDODGN = 1 53C ... always do diagonalization in this version 54 NROOT = MIN(IPRINT,NCONF,MXPRDT) 55C ... temporary analyze specification 56C 57C Core allocation: 58C 59 KHAMIL = 1 60 IF (IHSYM .EQ. 1) THEN 61 KFCAC2 = KHAMIL + (MXPRDT**2 + MXPRDT)/2 62 ELSE 63 KFCAC2 = KHAMIL + MXPRDT**2 64 END IF 65 KRJ = KFCAC2 + N2ASHX 66 KRK = KRJ + N2ASHX 67 KWHAM = KRK + N2ASHX 68 IF (IDODGN .EQ. 1) THEN 69 KEIGVL = KWHAM 70 KEIGVC = KEIGVL + MXPRDT 71 KWHAM = KEIGVC + MXPRDT ** 2 72 ELSE 73 KEIGVL = 999 999 999 74 KEIGVC = 999 999 999 75 END IF 76C 77 CALL DSPTSI(NASHT,FCAC,WRK(KFCAC2)) 78 CALL GETFIJ(WRK(KRJ),WRK(KRK),H2AC) 79 IF (MCTYPE .EQ. 2) THEN 80C ... reorder from Sirius to Lunar order for RAS 81 CALL REORMT(WRK(KFCAC2),WRK(KWHAM),NASHT,NASHT, 82 & XNDXCI(KLTSOB),XNDXCI(KLTSOB)) 83 CALL DCOPY(N2ASHX,WRK(KWHAM),1,WRK(KFCAC2),1) 84 CALL REORMT(WRK(KRJ),WRK(KWHAM),NASHT,NASHT, 85 & XNDXCI(KLTSOB),XNDXCI(KLTSOB)) 86 CALL DCOPY(N2ASHX,WRK(KWHAM),1,WRK(KRJ),1) 87 CALL REORMT(WRK(KRK),WRK(KWHAM),NASHT,NASHT, 88 & XNDXCI(KLTSOB),XNDXCI(KLTSOB)) 89 CALL DCOPY(N2ASHX,WRK(KWHAM),1,WRK(KRK),1) 90 END IF 91 CALL PHPDET(CIDIA,IWAY,MXPRDT,XNDXCI,WRK,KWHAM,LWRK, 92 & ICSYM,IHSYM,WRK(KHAMIL), 93 & WRK(KFCAC2),WRK(KRJ),WRK(KRK),NASHT,MULD2H, 94 & ECORE,ICOMBI,PSIGN,NPRDET,WRK(KEIGVC),WRK(KEIGVL), 95 & IDODGN,NCONF,IDET,NROOT,H2AC,IPRINT) 96C CALL PHPDET(CIDIA,IWAY,MXPRDT,XNDXCI,WORK,KFREE,LFREE, 97C & ICSYM,IHSYM,HAMIL,ONEBOD,RJ,RK,NORB,SYMPRO, 98C & ECORE,ICOMBI,PSIGN,NPRDET,EIGVC,EIGVL, 99C & IDODGN,NDET,IDET,NROOT,FIJKL,NTEST) 100C 101 IF (NPRDET .NE. MXPRDT) THEN 102 IF (IPRINT .GT. 0) THEN 103 WRITE (LUPRI,*) 'Number of explicit determinants' 104 WRITE (LUPRI,*) ' -- requested maximum',MXPRDT 105 WRITE (LUPRI,*) ' -- actually used ',NPRDET 106 END IF 107 MXPRDT = NPRDET 108 END IF 109 RETURN 110 END 111C /* Deck phpdet */ 112 SUBROUTINE PHPDET(CIDIA,IWAY,MXPRDT,XNDXCI,WORK,KFREE,LFREE, 113 & ICSYM,IHSYM,HAMIL,ONEBOD,RJ,RK,NORB,SYMPRO, 114 & ECORE,ICOMBI,PSIGN,NPRDET,EIGVC,EIGVL, 115 & IDODGN,NDET,IDET,NROOT,FIJKL,NTEST) 116C 117C.. SELECT A SUBSPACE AND CONSTRUCT HAMLITONIAN MATRIX IN THIS SUBSPACE 118C 119C PARAMETERS : 120C CIDIA : CI diagonal (input, DESTROYED if nroot.ne.0 ) 121C IWAY : flag for selecting subspace ( = 1 normal see below) 122C MXPRDT : largest allowed dimension of subspace 123C XNDXCI : string information 124C WORK : scratch space, length at least ?? 125C KFREE : first free element in work ( = 1 in normal cases ) 126C LFREE : dimension WORK(LFREE) 127C ICSYM : symmetry of CI space 128C IHSYM : = 0 : calculate complete subspace hamiltonian matrix 129C = 1 : only lower part of hamiltonian matrix calculated 130C HAMIL : space for hamilton matrix. DESTROYED if IDOGN .ne. 0 131C ONEBOD : onebody matrix in square form (FI for active orbitals ) 132C RJ : Coulomb integrals in square form 133C RK : Exchange integrals in square form 134C NORB : number of active orbitals 135C SYMPRO : 8*8 D2h table 136C ECORE : core energy 137C ICOMBI : = 0 : use determinant basis 138C = 1 : use combination basis 139C PSIGN : permutation sign for combinations ( only ICOMBI = 1 ) 140C NPRDET : Dimension of subspace used (OUTPUT ! ) 141C EIGVC : eigenvectors stored in column form ( IDOGN .NE. 0 ) 142C EIGVL : eigenvalues of subspace hamiltonian( IDOGN .NE. 0 ) 143C IDOGN : .ne. 0 : diagonalize subspace matrix 144C IDET : POINTER : IDET(I) gives adress of supspace element i in 145C FULL space ( output ). If IWAY = 3, also input. 146C NDET : dimension of full space 147C NROOT : number of roots analyzed with anaci ( IDOGN .ne. 0) 148C FIJKL : array relayed to H2TUVX routine in DIHDJ 149C NTEST : print level 150C 151C 152C.. SUBSPACE : 153C THE SUBSPACE IS CHOOSEN AS 154C IWAY = 1 : CHOOSE THE LOWEST VALUES OF 155C THE CI DIAGONAL.THE NUMBER OF DETERMINANTS IS 156C CHOOSEN SO NO DEGENERATE LEVELS ARE SPLITTED . 157C THE NUMBER OF DETERMINANTS USED , NPRDET, CAN 158C THEREFORE BE LOWER THAN MXPRDT. 159C IWAY = 2 : CHOOSE THE FIRST MXPRDT DETERMINANTS 160C STUPID, BUT CONVENIENT FOR TESTING 161C IWAY = 3 : Use input list in IDET(1:MXPRDT). List is not checked. 162C IWAY = 4 : choose the elements with highest absolut value 163C in CIDIA. (e.g. sigma vector for HF determinant). 164C 165C IDET CONTAINS ADRESSES OF ELEMENTS CHOSEN 166C 167C IDODGN .GT. 0 : DIAGONALIZE CONSTRUCTED HAMILTON MATRIX. 168C EIGVL CONTAINS EIGENVALUES ON RETURN 169C EIGVC CONTAINS EIGENVECTORS(COLUMNS) ON RETURN 170C 171C JEPPE OLSEN JANUARY 1989 172C 173#include "implicit.h" 174#include "priunit.h" 175 DIMENSION XNDXCI(*),WORK(*),HAMIL(*),CIDIA(*) 176 DIMENSION ONEBOD(*),RJ(*),RK(*),FIJKL(*) 177 DIMENSION IDET(MXPRDT), EIGVL(MXPRDT) , EIGVC(MXPRDT ** 2 ) 178 INTEGER SYMPRO(8,8) 179C 180C Used from common blocks: 181C STRNUM : ... 182C DETBAS : KLTSOB 183C 184#include "mxpdim.h" 185#include "strnum.h" 186#include "detbas.h" 187C 188 KFREEO = KFREE 189C 190 IF (NTEST .GT. 2) WRITE (LUPRI,*) '--- Output from PHPDET ---' 191 IF (NTEST .GT. 3) THEN 192 IF( ICOMBI .EQ. 0 ) THEN 193 WRITE(LUPRI,*) ' DETERMINANTS IN USE ' 194 ELSE 195 WRITE(LUPRI,*) ' COMBINATIONS IN USE ' 196 END IF 197 WRITE(LUPRI,*) ' MXPRDT ',MXPRDT 198 END IF 199C 200C.. 0 : SELECT SUBSPACE 201C 202 IF( IWAY .EQ. 1) THEN 203C. FIND NUMBER OF DETERMINANTS LESS OR EQUAL TOP MXPRDT 204C THAT DOES NOT SEPARATE DEGENERATE PAIRS . 205C FNDMN3(VEC,NDIM,MXELMN,IPLACE,NELMN,IPRT,THRES) 206 CALL FNDMN3(CIDIA,NDET,MXPRDT,IDET,NPRDET,NTEST,1.0D-3) 207 ELSE IF ( IWAY .EQ. 2 ) THEN 208 NPRDET = MIN(MXPRDT,NDET) 209 CALL ISTVC2(IDET,0,1,NPRDET) 210 ELSE IF ( IWAY .EQ. 3 ) THEN 211 NPRDET = MXPRDT 212C ... we assume caller supplies valid IDET(1:MXPRDT) and 213C MXPRDT .LE. NDET 214 ELSE IF ( IWAY .EQ. 4 ) THEN 215 CALL FNDAMX(CIDIA,NDET,MXPRDT,IDET,NPRDET,NTEST,1.0D-3) 216 ELSE 217 WRITE (LUPRI,*) 'Fatal error in PHPDET, illegal IWAY =',IWAY 218 CALL QUIT('Fatal error in PHPDET, illegal IWAY.') 219 END IF 220 IF (NTEST .GT. 10) WRITE(LUPRI,*) ' IWAY,NPRDET .... ', IWAY, 221 & NPRDET 222 IF (NTEST .GE. 10) THEN 223 WRITE(LUPRI,*) ' IDET IN PHPDET, IWAY = ',IWAY 224 CALL IWRTMA(IDET,1,NPRDET,1,NPRDET) 225 END IF 226C 227C.. 1 : LOCAL MEMORY 228C 229C SPACE FOR STRINGS OF SUBSPACE DETERMINANTS 230 CALL MEMADD(KIAST2,NPRDET*NAEL,KFREE,1) 231 CALL MEMADD(KIBST2,NPRDET*NBEL,KFREE,1) 232C FOR COMBINATIONS AN EXTRA SET IS NEEDED 233 IF( ICOMBI.NE. 0 ) THEN 234 CALL MEMADD(KIAST3,NPRDET*NAEL,KFREE,1) 235 CALL MEMADD(KIBST3,NPRDET*NBEL,KFREE,1) 236 ELSE 237 KIAST3 = KIAST2 238 KIBST3 = KIBST2 239 END IF 240C SCRATCH SPACE FOR DIHDJ 241 LSCR = 4 * NORB + NPRDET 242 CALL MEMADD(KSCR, LSCR,KFREE,1) 243 IF (KFREE .GT. LFREE) THEN 244 WRITE (LUPRI,*) 'Fatal error, insufficient memory in PHPDET' 245 WRITE (LUPRI,*) 'Need:',KFREE,' Available:',LFREE 246 CALL QUIT('Fatal error, insufficient memory in PHPDET') 247 END IF 248C 249C.. 3 : OBTAIN STRINGS CORRESPONDING TO SUBSPACE DETERMINANTS 250C 251 CALL STRFDT(NPRDET,IDET,WORK(KIAST2),WORK(KIBST2), 252 & ICSYM,SYMPRO,XNDXCI,ICOMBI) 253 IF( ICOMBI .NE. 0 ) THEN 254 CALL ICOPVE(WORK(KIAST2),WORK(KIAST3),NPRDET*NAEL) 255 CALL ICOPVE(WORK(KIBST2),WORK(KIBST3),NPRDET*NBEL) 256 END IF 257C 258C.. 4 : OBTAIN CORRESPONDING HAMILTONIAN 259C 260 CALL DIHDJ(WORK(KIAST2),WORK(KIBST2),NPRDET,WORK(KIAST3), 261 & WORK(KIBST3),NPRDET,NAEL,NBEL,WORK(KSCR),LSCR, 262 & NORB,ONEBOD,RJ,RK,FIJKL,XNDXCI(KLTSOB), 263 & HAMIL,IHSYM,ECORE,ICOMBI,PSIGN) 264C 265C.. 5 : DIAGONALIZE HAMIL ( ASSUMES IHSYM .NE. 0 USED ) 266C 267 IF( IDODGN .GT. 0 ) THEN 268 MROOT = MIN(NPRDET,NROOT) 269 CALL EIGEN(HAMIL,EIGVC,NPRDET,0,1) 270C. EXTRACT EIGENVALUES 271 CALL XTRCDI(HAMIL,EIGVL,NPRDET,1) 272 IF( NTEST.GE.7 ) THEN 273 WRITE(LUPRI,'(/A)') ' EIGENVALUES OF SUBSPACE HAMILTONIAN ' 274 CALL WRTMAT(EIGVL,1,NPRDET,1,NPRDET,0) 275 IF( NTEST .GE. 20 ) THEN 276 WRITE(LUPRI,'(/A)') ' EIGENVECTORS OF SUBSPACE HAMILTONIAN ' 277 CALL OUTPUT(EIGVC,1,NPRDET,1,NPRDET,NPRDET,NPRDET,1,LUPRI) 278 ELSE IF(NTEST .GE. 10 .AND. MROOT .GT. 0) THEN 279 WRITE(LUPRI,'(/A)') ' EIGENVECTORS OF SUBSPACE HAMILTONIAN ' 280 CALL OUTPUT(EIGVC,1,NPRDET,1,MROOT,NPRDET,NPRDET,1,LUPRI) 281 END IF 282 END IF 283 IF (NTEST.GT.2) WRITE(LUPRI,*) ' RANGE OF EIGENVALUES IN '// 284 & 'SUBSPACE ',EIGVL(1),EIGVL(NPRDET) 285C 286C .. 6 : ANALYZE MROOT APPROXIMATIVE EIGENVECTORS 287C 288 IF (NTEST .GE. 6) THEN 289 DO 1869 IROOT = 1, MROOT 290 IOFF = (IROOT-1)*NPRDET + 1 291 CALL SETVEC(CIDIA,0.0D0,NDET) 292 CALL SCAVEC(CIDIA,EIGVC(IOFF),IDET,NPRDET) 293 WRITE(LUPRI,'(/A,I3)') ' Information about root ... ',IROOT 294 WRITE(LUPRI,'(A)') ' ******************************' 295 WRITE(LUPRI,'(/A,1P,E15.8)') ' Energy .... ',EIGVL(IROOT) 296 CUTOFF = 0.1D0/SQRT(10.0D0) 297 CALL ANACIN(CIDIA,ICSYM,CUTOFF,100,XNDXCI,SYMPRO,6,0) 298C CALL ANACIN(CIVEC,ICSYM,THRES,MAXTRM,XNDXCI,SYMPRO,IOUT,INCSFB) 299 1869 CONTINUE 300 END IF 301 END IF 302C 303 KFREE = KFREEO 304 RETURN 305 END 306C /* Deck strfdt */ 307 SUBROUTINE STRFDT(NDET,IDET,IASTR2,IBSTR2,ICSYM,SYMPRO, 308 & XNDXCI,JCOMBI) 309C 310C NDET DETERMINANTS ARE GIVEN IN ARRAY IDET 311C OBTAIN THE CORRESPONDING ALPHA AND BETASTRINGS AND STORE THEM 312C IN IASTR2 AND IBSTR2 313C 314C NOTE ICOMBI IS INTRODUCED THROUGH CIINFO IN THIS IMPLEMENTATION 315C 316#include "implicit.h" 317 INTEGER SYMPRO(8,8) 318 DIMENSION IDET(*),IASTR2(*),IBSTR2(*),XNDXCI(*) 319C 320#include "mxpdim.h" 321#include "ciinfo.h" 322#include "detbas.h" 323#include "strnum.h" 324C 325 CALL STFDT2(NDET,IDET,IASTR2,IBSTR2,MAXSYM,NOCTPA,NOCTPB, 326 & XNDXCI(KNSSOA),XNDXCI(KNSSOB),XNDXCI(KISSOA), 327 & XNDXCI(KISSOB),XNDXCI(KIOCOC),NAEL,NBEL, 328 & SYMPRO,XNDXCI(KIASTR),XNDXCI(KIBSTR),ICSYM,ICOMBI) 329C SUBROUTINE STFDT2(NDET,IDET,IASTR2,IBSTR2,MAXSYM,NOCTPA,NOCTPB, 330C & NBEL,SYMPRO,IASTR,IBSTR,ISYM,ICOMBI) 331C 332 RETURN 333 END 334C /* Deck dihdj */ 335 SUBROUTINE DIHDJ(IASTR,IBSTR,NIDET,JASTR,JBSTR,NJDET,NAEL,NBEL, 336 & IWORK,LWORK,NORB,ONEBOD,RJ,RK,FIJKL,ILTSOB, 337 & HAMIL,ISYM,ECORE,ICOMBI,PSIGN) 338C 339C A SET OF DETERMINANTS IA DEFINED BY ALPHA AND BETASTRINGS 340C IASTR,IBSTR AND ANOTHER SET OF DETERMINATS DEFINED BY STRINGS 341C JASTR AND JBSTR ARE GIVEN . OBTAIN CORRESPONDING HAMILTONIAN MATRIX 342C 343C IF ICOMBI .NE. 0 COMBINATIONS ARE USED FOR ALPHA AND BETA STRING 344C THAT DIFFERS : 345C 1/SQRT(2) * ( |I1A I2B| + PSIGN * |I2A I1B| ) 346C 347C IF ISYM .EQ. 0 FULL HAMILTONIAN IS CONSTRUCTED 348C IF ISYM .NE. 0 LOWER HALF OF HAMILTONIAN IS CONSTRUCTED 349C 350C JEPPE OLSEN JANUARY 1989 351C 352#include "implicit.h" 353#include "priunit.h" 354 DIMENSION IASTR(NAEL,*),IBSTR(NBEL,*) 355 DIMENSION JASTR(NAEL,*),JBSTR(NBEL,*) 356C 357 DIMENSION IWORK(*), HAMIL(*), ONEBOD(NORB,NORB), 358 * RJ(NORB,NORB), RK(NORB,NORB), FIJKL(*), ILTSOB(*) 359C 360C SCRATCH SPACE 361C 362C .. 1 : EXPANSION OF ALPHA AND BETA STRINGS OF TYPE I 363C 364C? WRITE(LUPRI,*) ' MISTER DIHDJ SPEAKING ICOMBI = ', ICOMBI 365 KLFREE = 1 366 KLIAE = KLFREE 367 KLFREE = KLIAE + NORB 368 KLIBE = KLFREE 369 KLFREE = KLIBE + NORB 370C 371 KLJAE = KLFREE 372 KLFREE = KLJAE + NORB 373 KLJBE = KLFREE 374 KLFREE = KLJBE + NORB 375 IF( ICOMBI .NE. 0 ) THEN 376 KLIAB = KLFREE 377 KLFREE = KLFREE + NIDET 378 END IF 379C? WRITE (LUPRI,*) ' PSIGN ', PSIGN 380C 381 IF( ICOMBI .NE. 0 ) THEN 382C SET UP ARRAY COMBARING ALPHA AND BETA STRINGS IN IDET LIST 383 DO 13 IDET = 1, NIDET 384 IAEQIB = 1 385 DO 14 IEL = 1, NAEL 386 IF(IASTR(IEL,IDET) .NE. IBSTR(IEL,IDET))IAEQIB = 0 387 14 CONTINUE 388 IWORK(KLIAB-1+IDET) = IAEQIB 389 13 CONTINUE 390C? WRITE(LUPRI,*) ' IAEQIB ARRAY : ' 391C? CALL IWRTMA(IWORK(KLIAB),1,NIDET,1,NIDET) 392 END IF 393C 394 IF( ISYM .EQ. 0 ) THEN 395 LHAMIL = NIDET*NJDET 396 ELSE 397 LHAMIL = NIDET*(NIDET+1) / 2 398 END IF 399 CALL SETVEC(HAMIL,0.0D0,LHAMIL) 400C 401 NTERMS= 0 402 NDIF0 = 0 403 NDIF1 = 0 404 NDIF2 = 0 405C.. LOOP OVER J DETERMINANTS 406C 407 NTEST = 0 408 DO 1000 JDET = 1,NJDET 409 IF( NTEST .GE. 20 ) WRITE(LUPRI,*) ' ****** 1000 JDET ', JDET 410C 411C EXPAND JDET 412 CALL ISETVC(IWORK(KLJAE),0,NORB) 413 CALL ISETVC(IWORK(KLJBE),0,NORB) 414C 415 IF( ICOMBI .NE. 0 ) THEN 416 JAEQJB = 1 417 DO 32 IEL = 1, NAEL 418 IF(JASTR(IEL,JDET) .NE. JBSTR(IEL,JDET))JAEQJB = 0 419 32 CONTINUE 420C? WRITE(LUPRI,*) ' JAEQJB ', JAEQJB 421 END IF 422C 423 DO 40 IAEL = 1, NAEL 424 IWORK(KLJAE-1+JASTR(IAEL,JDET) ) = 1 425 40 CONTINUE 426C 427 DO 50 IBEL = 1, NBEL 428 IWORK(KLJBE-1+JBSTR(IBEL,JDET) ) = 1 429 50 CONTINUE 430C 431 IF( NTEST .GE. 10 ) THEN 432 WRITE(LUPRI,*) ' LOOP 1000 JDET = ',JDET 433C WRITE(LUPRI,*) ' JASTR AND JBSTR ' 434C CALL IWRTMA(JASTR(1,JDET),1,NAEL,1,NAEL) 435C CALL IWRTMA(JBSTR(1,JDET),1,NBEL,1,NBEL) 436C WRITE(LUPRI,*) ' EXPANDED ALPHA AND BETA STRING ' 437C CALL IWRTMA(IWORK(KLJAE),1,NORB,1,NORB) 438C CALL IWRTMA(IWORK(KLJBE),1,NORB,1,NORB) 439 END IF 440C 441 IF( ISYM .EQ. 0 ) THEN 442 MINI = 1 443 ELSE 444 MINI = JDET 445 END IF 446 DO 900 IDET = MINI, NIDET 447C? WRITE(LUPRI,*) ' LOOP 900 IDET .... ',IDET 448 IF( ICOMBI .EQ. 0 ) THEN 449 NLOOP = 1 450 ELSE 451 IAEQIB = IWORK(KLIAB-1+IDET) 452 IF(IAEQIB+JAEQJB .EQ. 0 ) THEN 453 NLOOP = 2 454 ELSE 455 NLOOP = 1 456 END IF 457 END IF 458C 459 DO 899 ILOOP = 1, NLOOP 460 NTERMS = NTERMS + 1 461C? WRITE(LUPRI,*) ' 899 : ILOOP ' , ILOOP 462C 463C.. COMPARE DETERMINANTS 464C 465C SWAP IA AND IB FOR SECOND PART OF COMBINATIONS 466 IF( ILOOP .EQ. 2 ) 467 & CALL ISWPVE(IASTR(1,IDET),IBSTR(1,IDET),NAEL) 468C 469 NACM = 0 470 DO 61 IAEL = 1, NAEL 471 NACM = NACM + IWORK(KLJAE-1+IASTR(IAEL,IDET)) 472 61 CONTINUE 473C 474 NBCM = 0 475 DO 62 IBEL = 1, NBEL 476 NBCM = NBCM + IWORK(KLJBE-1+IBSTR(IBEL,IDET)) 477 62 CONTINUE 478C 479 NADIF = NAEL-NACM 480 NBDIF = NBEL-NBCM 481 IF( NTEST .GE. 10 ) THEN 482 WRITE(LUPRI,*) ' LOOP 900 IDET ',IDET 483 WRITE(LUPRI,*) ' COMPARISON , NADIF , NBDIF ', NADIF,NBDIF 484 END IF 485C 486 IF(NADIF+NBDIF .GT. 2 ) GOTO 898 487C 488C 489C FACTOR FOR COMBINATIONS 490 IF( ICOMBI .EQ. 0 ) THEN 491 CONST = 1.0D0 492 ELSE 493 IF((JAEQJB +IAEQIB) .EQ.2 ) THEN 494 CONST = 1.0D0 495 ELSE IF( (JAEQJB+IAEQIB) .EQ. 1 ) THEN 496COLD CONST = 1.0D0/SQRT(2.0D0)*(1.0D0+PSIGN) 497 IF (PSIGN .EQ. -1.0D0) GO TO 898 498C ... no contribution between singlet and triplet. 499 CONST = SQRT(2.0D0) 500 ELSE IF( (JAEQJB+IAEQIB) .EQ. 0 ) THEN 501 IF( ILOOP .EQ. 1) THEN 502 CONST = 1.0D0 503 ELSE 504 CONST = PSIGN 505 END IF 506 END IF 507 END IF 508C 509C.. FIND DIFFERING ORBITALS AND SIGN FOR PERMUTATION 510C 511C EXPAND IDET 512 CALL ISETVC(IWORK(KLIAE),0,NORB) 513 CALL ISETVC(IWORK(KLIBE),0,NORB) 514C 515 DO 42 IAEL = 1, NAEL 516 IWORK(KLIAE-1+IASTR(IAEL,IDET) ) = 1 517 42 CONTINUE 518C 519 DO 52 IBEL = 1, NBEL 520 IWORK(KLIBE-1+IBSTR(IBEL,IDET) ) = 1 521 52 CONTINUE 522C 523 IF(NADIF .EQ. 1 ) THEN 524 DO 120 IAEL = 1,NAEL 525 IF(IWORK(KLJAE-1+IASTR(IAEL,IDET)).EQ.0) THEN 526 IA = IASTR(IAEL,IDET) 527 IEL1 = IAEL 528 GOTO 121 529 END IF 530 120 CONTINUE 531 121 CONTINUE 532C 533 DO 130 JAEL = 1,NAEL 534 IF(IWORK(KLIAE-1+JASTR(JAEL,JDET)).EQ.0) THEN 535 JA = JASTR(JAEL,JDET) 536 JEL1 = JAEL 537 GOTO 131 538 END IF 539 130 CONTINUE 540 131 CONTINUE 541 SIGNA = (-1)**(JEL1+IEL1) 542C? WRITE(LUPRI,*) ' IA JA SIGNA... ',IA,JA,SIGNA 543C 544 END IF 545C 546 IF(NBDIF .EQ. 1 ) THEN 547 DO 220 IBEL = 1,NBEL 548 IF(IWORK(KLJBE-1+IBSTR(IBEL,IDET)).EQ.0) THEN 549 IB = IBSTR(IBEL,IDET) 550 IEL1 = IBEL 551 GOTO 221 552 END IF 553 220 CONTINUE 554 221 CONTINUE 555C 556 DO 230 JBEL = 1,NBEL 557 IF(IWORK(KLIBE-1+JBSTR(JBEL,JDET)).EQ.0) THEN 558 JB = JBSTR(JBEL,JDET) 559 JEL1 = JBEL 560 GOTO 231 561 END IF 562 230 CONTINUE 563 231 CONTINUE 564 SIGNB = (-1)**(JEL1+IEL1) 565C? WRITE(LUPRI,*) ' IB JB SIGNB... ',IB,JB,SIGNB 566C 567 END IF 568 IF(NADIF .EQ. 2 ) THEN 569 IDIFF = 0 570 DO 320 IAEL = 1,NAEL 571 IF(IWORK(KLJAE-1+IASTR(IAEL,IDET)).EQ.0) THEN 572 IF( IDIFF .EQ. 0 ) THEN 573 IDIFF = 1 574 I1 = IASTR(IAEL,IDET) 575 IPERM = IAEL 576 ELSE 577 I2 = IASTR(IAEL,IDET) 578 IPERM = IAEL + IPERM 579 GOTO 321 580 END IF 581 END IF 582 320 CONTINUE 583 321 CONTINUE 584C 585 JDIFF = 0 586 DO 330 JAEL = 1,NAEL 587 IF(IWORK(KLIAE-1+JASTR(JAEL,JDET)).EQ.0) THEN 588 IF( JDIFF .EQ. 0 ) THEN 589 JDIFF = 1 590 J1 = JASTR(JAEL,JDET) 591 JPERM = JAEL 592 ELSE 593 J2 = JASTR(JAEL,JDET) 594 JPERM = JAEL + JPERM 595 GOTO 331 596 END IF 597 END IF 598 330 CONTINUE 599 331 CONTINUE 600 SIGN = (-1)**(IPERM+JPERM) 601C 602 END IF 603C 604 IF(NBDIF .EQ. 2 ) THEN 605 IDIFF = 0 606 DO 420 IBEL = 1,NBEL 607 IF(IWORK(KLJBE-1+IBSTR(IBEL,IDET)).EQ.0) THEN 608 IF( IDIFF .EQ. 0 ) THEN 609 IDIFF = 1 610 I1 = IBSTR(IBEL,IDET) 611 IPERM = IBEL 612 ELSE 613 I2 = IBSTR(IBEL,IDET) 614 IPERM = IBEL + IPERM 615 GOTO 421 616 END IF 617 END IF 618 420 CONTINUE 619 421 CONTINUE 620C 621 JDIFF = 0 622 DO 430 JBEL = 1,NBEL 623 IF(IWORK(KLIBE-1+JBSTR(JBEL,JDET)).EQ.0) THEN 624 IF( JDIFF .EQ. 0 ) THEN 625 JDIFF = 1 626 J1 = JBSTR(JBEL,JDET) 627 JPERM = JBEL 628 ELSE 629 J2 = JBSTR(JBEL,JDET) 630 JPERM = JBEL + JPERM 631 GOTO 431 632 END IF 633 END IF 634 430 CONTINUE 635 431 CONTINUE 636 SIGN = (-1)**(IPERM+JPERM) 637C 638 END IF 639C 640C OBTAIN VALUE OF HAMILTONIAN ELEMENT 641C 642 IF( NADIF .EQ. 2 .OR. NBDIF .EQ. 2 ) THEN 643 NDIF2 = NDIF2 + 1 644C SIGN * (I1 J1 | I2 J2 ) - ( I1 J2 | I2 J1 ) 645 XVAL = SIGN*( H2TUVX(I1,J1,I2,J2,FIJKL,ILTSOB) 646 * -H2TUVX(I1,J2,I2,J1,FIJKL,ILTSOB) ) 647 ELSE IF( NADIF .EQ. 1 .AND. NBDIF .EQ. 1 ) THEN 648 NDIF2 = NDIF2 + 1 649C SIGN * (IA JA | IB JB ) 650C? WRITE(LUPRI,*) ' IA IB JA JB ', IA,IB,JA,JB 651 XVAL = SIGNA*SIGNB* H2TUVX(IA,JA,IB,JB,FIJKL,ILTSOB) 652C? WRITE(LUPRI,*) ' SIGNA SIGNB XVAL ',SIGNA,SIGNB,XVAL 653 ELSE IF( NADIF .EQ. 1 .AND. NBDIF .EQ. 0 .OR. 654 & NADIF .EQ. 0 .AND. NBDIF .EQ. 1 )THEN 655 NDIF1 = NDIF1 + 1 656C SIGN * 657C( H(I1 J1 ) + 658C (SUM OVER ORBITALS OF BOTH SPIN TYPES ( I1 J1 | JORB JORB ) 659C -(SUM OVER ORBITALS OF DIFFERING SPIN TYPE ( I1 JORB | JORB J1 ) ) 660 IF( NADIF .EQ. 1 ) THEN 661 I1 = IA 662 J1 = JA 663 SIGN = SIGNA 664 ELSE 665 I1 = IB 666 J1 = JB 667 SIGN = SIGNB 668 END IF 669C? WRITE(LUPRI,*) ' ONE DIFF I1 J1 SIGN : ',I1,J1,SIGN 670C 671 XVAL = ONEBOD(I1,J1) 672 DO 520 JAEL = 1, NAEL 673 JORB = JASTR(JAEL,JDET) 674 XVAL = XVAL + H2TUVX(I1,J1,JORB,JORB,FIJKL,ILTSOB) 675 520 CONTINUE 676 DO 521 JBEL = 1, NBEL 677 JORB = JBSTR(JBEL,JDET) 678 XVAL = XVAL + H2TUVX(I1,J1,JORB,JORB,FIJKL,ILTSOB) 679 521 CONTINUE 680 IF( NADIF .EQ. 1 ) THEN 681 DO 522 JAEL = 1, NAEL 682 JORB = JASTR(JAEL,JDET) 683 XVAL = XVAL - H2TUVX(I1,JORB,JORB,J1,FIJKL,ILTSOB) 684 522 CONTINUE 685 ELSE 686 DO 523 JBEL = 1, NBEL 687 JORB = JBSTR(JBEL,JDET) 688 XVAL = XVAL - H2TUVX(I1,JORB,JORB,J1,FIJKL,ILTSOB) 689 523 CONTINUE 690 END IF 691 XVAL = XVAL * SIGN 692 ELSE IF( NADIF .EQ. 0 .AND. NBDIF .EQ. 0 ) THEN 693 NDIF0 = NDIF0 + 1 694C SUM(I,J OF JDET) H(I,J) + (I I | J J ) - (I J | J I ) 695C 696 XVAL = ECORE 697 DO 650 IAB = 1, 2 698 IF(IAB .EQ. 1 ) THEN 699 NIABEL = NAEL 700 ELSE 701 NIABEL = NBEL 702 END IF 703 DO 640 JAB = 1, 2 704 IF(JAB .EQ. 1 ) THEN 705 NJABEL = NAEL 706 ELSE 707 NJABEL = NBEL 708 END IF 709 DO 630 IEL = 1, NIABEL 710 IF( IAB .EQ. 1 ) THEN 711 IORB = IASTR(IEL,IDET) 712 ELSE 713 IORB = IBSTR(IEL,IDET) 714 END IF 715 IF(IAB .EQ. JAB ) XVAL = XVAL + ONEBOD(IORB,IORB) 716 IORB = IORB 717 DO 620 JEL = 1, NJABEL 718 IF( JAB .EQ. 1 ) THEN 719 JORB = IASTR(JEL,IDET) 720 ELSE 721 JORB = IBSTR(JEL,IDET) 722 END IF 723 XVAL = XVAL + 0.5D0*RJ(IORB,JORB) 724 IF( IAB . EQ. JAB ) 725 & XVAL = XVAL - 0.5D0*RK(IORB,JORB) 726 620 CONTINUE 727 630 CONTINUE 728 640 CONTINUE 729 650 CONTINUE 730 END IF 731C 732C? WRITE(LUPRI,*) ' CONST XVAL ', CONST ,XVAL 733 IF( ISYM .EQ. 0 ) THEN 734 HAMIL((JDET-1)*NIDET+IDET) = 735 & HAMIL((JDET-1)*NIDET+IDET) + CONST * XVAL 736 ELSE 737 HAMIL((IDET-1)*IDET/2 + JDET ) = 738 & HAMIL((IDET-1)*IDET/2 + JDET ) + CONST * XVAL 739 END IF 740C RESTORE ORDER !!! 741 898 IF( ILOOP .EQ. 2 ) 742 & CALL ISWPVE(IASTR(1,IDET),IBSTR(1,IDET),NAEL) 743 899 CONTINUE 744 900 CONTINUE 745 1000 CONTINUE 746C 747 IF (NTEST .GE.2) THEN 748 WRITE(LUPRI,*) 749 WRITE(LUPRI,*) 750 &' Number of elements differing by 0 excitation.. ',NDIF0 751C 752 WRITE(LUPRI,*) 753 &' Number of elements differing by 1 excitation.. ',NDIF1 754C 755 WRITE(LUPRI,*) 756 &' Number of elements differing by 2 excitation.. ',NDIF2 757C 758 WRITE(LUPRI,*) 759 &' Number of vanishing elments ', 760 & NTERMS - NDIF0 - NDIF1 - NDIF2 761 END IF 762 IF( NTEST .GT. 2 ) THEN 763 WRITE(LUPRI,'(/A)') ' HAMILTONIAN MATRIX ' 764 IF( ISYM .EQ. 0 ) THEN 765 CALL WRTMAT(HAMIL,NIDET,NJDET,NIDET,NJDET,0) 766 ELSE 767 CALL PRSYM(HAMIL,NIDET) 768 END IF 769 END IF 770C 771C! CALL QUIT(' ENFORCED STOP AT END OF DIHDJ ') 772 RETURN 773 END 774C /* Deck stfdt2 */ 775 SUBROUTINE STFDT2(NDET,IDET,IASTR2,IBSTR2,MAXSYM,NOCTPA,NOCTPB, 776 & NSSOA,NSSOB,ISSOA,ISSOB,IOCOC,NAEL, 777 & NBEL,SYMPRO,IASTR,IBSTR,ISYM,ICOMBI) 778C 779C EXTRACT ALPHA AND BETASTRINGS FOR NDET DETERMINANTS GIVEN IN 780C IDET 781C 782C RAS VERSION 783C 784C JEPPE OLSEN , JANUARY 1989 785C 786#include "implicit.h" 787#include "priunit.h" 788 DIMENSION IDET(NDET),IASTR2(NAEL,NDET),IBSTR2(NBEL,NDET) 789C 790 DIMENSION NSSOA(NOCTPA,MAXSYM),NSSOB(NOCTPB,MAXSYM) 791 DIMENSION ISSOA(NOCTPA,MAXSYM),ISSOB(NOCTPB,MAXSYM) 792 DIMENSION IOCOC(NOCTPA,NOCTPB) 793 INTEGER SYMPRO(8,8) 794 DIMENSION IASTR(NAEL,*),IBSTR(NBEL,*) 795C 796 NFOUND = 0 797 LDET = 0 798 DO 1000 IASYM = 1, MAXSYM 799 IBSYM = SYMPRO(ISYM,IASYM) 800 DO 999 IATP = 1,NOCTPA 801 NIA = NSSOA(IATP,IASYM) 802 IF (NIA.EQ.0) GO TO 999 803 IAOFF = ISSOA(IATP,IASYM) 804 DO 901 IBTP = 1, NOCTPB 805 IF(IOCOC(IATP,IBTP).LE.0) GOTO 901 806 NIB = NSSOB(IBTP,IBSYM) 807 IBOFF = ISSOB(IBTP,IBSYM) 808 DO 900 IB = 1, NIB 809 DO 800 IA = 1, NIA 810 IBEFF = IBOFF-1+IB 811 IAEFF = IAOFF-1+IA 812C? WRITE(LUPRI,*) ' IAEFF IBEFF ' , IAEFF,IBEFF 813C THE FOLLOWING IS NEW LOGICAL ORDERING ! 814 IF(ICOMBI.NE.0..AND.IAEFF.LT.IBEFF) 815 & GOTO 800 816 LDET = LDET + 1 817 NEXT = 0 818 DO 500 I =1, NDET 819 500 IF(IDET(I) .EQ. LDET ) NEXT = I 820 IF(NEXT.NE.0 ) THEN 821 CALL ICOPVE(IASTR(1,IAEFF),IASTR2(1,NEXT),NAEL) 822 CALL ICOPVE(IBSTR(1,IBEFF),IBSTR2(1,NEXT),NBEL) 823 NFOUND = NFOUND + 1 824 IF(NFOUND .EQ. NDET ) GOTO 1001 825 END IF 826 800 CONTINUE 827 900 CONTINUE 828 901 CONTINUE 829 999 CONTINUE 830 1000 CONTINUE 831 1001 CONTINUE 832C 833 NTEST = 0 834 IF( NTEST .GE. 2 ) THEN 835 IF( ICOMBI .EQ. 0 ) THEN 836 WRITE(LUPRI,*) ' DETERMINANTS SELECTED ' 837 ELSE 838 WRITE(LUPRI,*) ' DETERMINANT COMBINATIONS SELECTED ' 839 END IF 840C 841 CALL IWRTMA(IDET,1,NDET,1,NDET) 842 WRITE(LUPRI,*) ' SELECTED ALPHA AND BETA STRINGS FROM STRDT1 ' 843 CALL IWRTMA(IASTR2,NAEL,NDET,NAEL,NDET) 844 WRITE(LUPRI,*) 845 CALL IWRTMA(IBSTR2,NBEL,NDET,NBEL,NDET) 846 END IF 847C 848 RETURN 849 END 850C /* Deck fndmn3 */ 851 SUBROUTINE FNDMN3(VEC,NDIM,MXELMN,IPLACE,NELMN,IPRT,THRES) 852C 853C Jeppe Olsen. Last revision 890205/900405 hjaaj. 854C 855C FIND MXELMN/NELMN LOWEST ELEMENTS IN VEC . 856C NELMN IS THE LARGEST NUMBER LOWER THAN MXELMN THAT DOES NOT 857C SPLIT DEGENERATE PAIRS 858C ORIGINAL PLACES OF THE LOWEST ELEMENTS ARE STORED IN IPLACE 859C 860#include "implicit.h" 861#include "priunit.h" 862 DIMENSION VEC(NDIM),IPLACE(*) 863C 864 PARAMETER ( D1 = 1.0D0 ) 865C 866C. FIRST OCCURANCE OF LOWEST ELEMENT AND LARGEST ELEMENT 867 XMIN = VEC(1) 868 XMAX = VEC(1) 869 IMIN = 1 870 IMAX = 1 871 DO 100 I = 2,NDIM 872 IF( VEC(I) .GT. XMAX) THEN 873 XMAX = VEC(I) 874 IMAX = I 875 ELSE IF( VEC(I) .LT. XMIN) THEN 876 XMIN = VEC(I) 877 IMIN = I 878 END IF 879 100 CONTINUE 880C 881 IF (IPRT.GT.0) WRITE(LUPRI,*) ' -- output from FNDMN3 --' 882 IF(IPRT .GE. 5 ) THEN 883 WRITE(LUPRI,*) ' LOWEST VALUE AND PLACE ',XMIN,IMIN 884 WRITE(LUPRI,*) ' HIGHST VALUE AND PLACE ',XMAX,IMAX 885 END IF 886C 887 IPLACE(1) = IMIN 888 NDEG = 1 889 IF(IPRT .GE. 5 ) WRITE(LUPRI,'(A,I8,1P,D15.5,I8,I4)') 890 & 'IELMNT XMIN IMIN NDEG ', 1,XMIN,IMIN,NDEG 891C 892 ITOP = MIN(NDIM,MXELMN+1) 893 DO 200 IELMNT = 2,ITOP 894 XMINPR = XMIN 895 IMINPR = IMIN 896 XMIN = XMAX + D1 897 IMIN = -1 898 DO 150 I = 1,NDIM 899 IF (VEC(I) .LT. XMIN) THEN 900 IF (VEC(I) .GT. XMINPR) THEN 901 IMIN = I 902 XMIN = VEC(I) 903 ELSE IF (VEC(I) .EQ. XMINPR .AND. I .GT. IMINPR) THEN 904 IMIN = I 905 XMIN = VEC(I) 906 GO TO 151 907 END IF 908 END IF 909 150 CONTINUE 910 151 CONTINUE 911 IF(XMIN-XMINPR .LE. THRES )THEN 912 NDEG = NDEG + 1 913 ELSE 914 NDEG = 1 915 END IF 916C 917 IF( IELMNT .LE. MXELMN ) IPLACE(IELMNT) = IMIN 918 IF( IPRT .GE. 5 ) WRITE(LUPRI,'(A,I8,1P,D15.5,I8,I4)') 919 & 'IELMNT XMIN IMIN NDEG ', IELMNT,XMIN,IMIN,NDEG 920 200 CONTINUE 921C 922C CHECK DEGENERACY ON LAST VALUE 923 IF(MXELMN .LT. NDIM ) THEN 924 NELMN = MXELMN+1-NDEG 925 ELSE 926 NELMN = NDIM 927 END IF 928 IF (IPRT .GT. 0) 929 & WRITE(LUPRI,*) ' NUMBER OF ELEMENTS OBTAINED IN FNDMN3 ',NELMN 930C 931C 932 IF( IPRT .GE. 3 ) THEN 933 WRITE(LUPRI,*) ' FROM FNDMN3 : ' 934 WRITE(LUPRI,*) ' PLACES OF LOWEST ELEMENTS ' 935 CALL IWRTMA(IPLACE,1,NELMN ,1,NELMN ) 936 END IF 937C 938 IF( IPRT .GE. 1 ) THEN 939 WRITE(LUPRI,*) 940 & ' MIN AND MAX IN SELECTED SUPSPACE ', 941 & VEC(IPLACE(1)),VEC(IPLACE(NELMN)) 942 END IF 943C 944 RETURN 945 END 946C /* Deck fndamx */ 947 SUBROUTINE FNDAMX(VEC,NDIM,MXELMN,IPLACE,NELMN,IPRT,THRES) 948C 949C 890205 hjaaj. Based on FNDMN3 by Jeppe Olsen. 950C Last revision 900405 hjaaj. 951C 952C FIND MXELMN/NELMN ELEMENTS IN VEC with highest absolute value. 953C NELMN IS THE LARGEST NUMBER LOWER THAN MXELMN THAT DOES NOT 954C SPLIT DEGENERATE PAIRS. 955C ORIGINAL PLACES OF THE LOWEST ELEMENTS ARE STORED IN IPLACE. 956C 957#include "implicit.h" 958#include "priunit.h" 959 DIMENSION VEC(NDIM),IPLACE(*) 960C 961 PARAMETER ( D1 = 1.0D0 ) 962C 963C. FIRST OCCURANCE OF LOWEST ELEMENT AND LARGEST ELEMENT 964 XMIN = ABS(VEC(1)) 965 XMAX = ABS(VEC(1)) 966 IMIN = 1 967 IMAX = 1 968 DO 100 I = 2,NDIM 969 AVECI = ABS(VEC(I)) 970 IF( AVECI .GT. XMAX) THEN 971 XMAX = AVECI 972 IMAX = I 973 ELSE IF( AVECI .LT. XMIN) THEN 974 XMIN = AVECI 975 IMIN = I 976 END IF 977 100 CONTINUE 978C 979 IF (IPRT.GT.0) WRITE(LUPRI,*) ' -- output from FNDAMX --' 980 IF(IPRT .GE. 5 ) THEN 981 WRITE(LUPRI,*) ' LOWEST VALUE AND PLACE ',XMIN,IMIN 982 WRITE(LUPRI,*) ' HIGHST VALUE AND PLACE ',XMAX,IMAX 983 END IF 984C 985 IPLACE(1) = IMAX 986 NDEG = 1 987 IF(IPRT .GE. 5 ) WRITE(LUPRI,'(A,I8,1P,D15.5,I8,I4)') 988 & 'IELMNT value IMAX NDEG ', 1,VEC(IMAX),IMAX,NDEG 989C 990 ITOP = MIN(NDIM,MXELMN+1) 991 DO 200 IELMNT = 2,ITOP 992 XMAXPR = XMAX 993 IMAXPR = IMAX 994 XMAX = -D1 995 IMAX = -1 996 DO 150 I = 1,NDIM 997 AVECI = ABS(VEC(I)) 998 IF (AVECI .GT. XMAX) THEN 999 IF (AVECI .LT. XMAXPR) THEN 1000 IMAX = I 1001 XMAX = AVECI 1002 ELSE IF (AVECI .EQ. XMAXPR .AND. I .GT. IMAXPR) THEN 1003 IMAX = I 1004 XMAX = AVECI 1005 GO TO 151 1006 END IF 1007 END IF 1008 150 CONTINUE 1009 151 CONTINUE 1010 IF(XMAXPR-XMAX .LE. THRES )THEN 1011 NDEG = NDEG + 1 1012 ELSE 1013 NDEG = 1 1014 END IF 1015C 1016 IF( IELMNT .LE. MXELMN ) IPLACE(IELMNT) = IMAX 1017 IF( IPRT .GE. 5 ) WRITE(LUPRI,'(A,I8,1P,D15.5,I8,I4)') 1018 & ' IELMNT value IMAX NDEG ',IELMNT,VEC(IMAX),IMAX,NDEG 1019 200 CONTINUE 1020C 1021C CHECK DEGENERACY ON LAST VALUE 1022 IF(MXELMN .LT. NDIM ) THEN 1023 NELMN = MXELMN+1-NDEG 1024 ELSE 1025 NELMN = NDIM 1026 END IF 1027 IF (IPRT .GT. 0) 1028 & WRITE(LUPRI,*) ' NUMBER OF ELEMENTS OBTAINED IN FNDAMX ',NELMN 1029C 1030C 1031 IF( IPRT .GE. 3 ) THEN 1032 WRITE(LUPRI,*) ' FROM FNDAMX : ' 1033 WRITE(LUPRI,*) ' PLACES OF ELEMENTS of highest absolut value' 1034 CALL IWRTMA(IPLACE,1,NELMN ,1,NELMN ) 1035 END IF 1036C 1037 IF( IPRT .GE. 1 ) THEN 1038 WRITE(LUPRI,*) 1039 & ' MAX AND MIN IN SELECTED SUPSPACE ', 1040 & VEC(IPLACE(1)),VEC(IPLACE(NELMN)) 1041 END IF 1042C 1043 RETURN 1044 END 1045C /* Deck tripk2 */ 1046 SUBROUTINE TRIPK2(AUTPAK,APAK,IWAY,MATDIM,NDIM,SIGN) 1047C 1048C 1049C.. REFORMATING BETWEEN LOWER TRIANGULAR PACKING 1050C AND FULL MATRIX FORM FOR A SYMMETRIC OR ANTI SYMMETRIC MATRIX 1051C 1052C IWAY = 1 : FULL TO PACKED 1053C LOWER HALF OF AUTPAK IS STORED IN APAK 1054C IWAY = 2 : PACKED TO FULL FORM 1055C APAK STORED IN LOWER HALF 1056C SIGN * APAK TRANSPOSED IS STORED IN UPPPER PART 1057C.. NOTE : COLUMN WISE STORAGE SCHEME IS USED FOR PACKED BLOCKS 1058#include "implicit.h" 1059#include "priunit.h" 1060 DIMENSION AUTPAK(MATDIM,MATDIM),APAK(*) 1061C 1062 IF( IWAY .EQ. 1 ) THEN 1063 IJ = 0 1064 DO 100 J = 1,NDIM 1065 DO 50 I = J , NDIM 1066 APAK(IJ+I) = AUTPAK(I,J) 1067 50 CONTINUE 1068 IJ = IJ +NDIM-J 1069 100 CONTINUE 1070 END IF 1071C 1072 IF( IWAY .EQ. 2 ) THEN 1073 IJ = 0 1074 DO 200 J = 1,NDIM 1075 DO 150 I = J,NDIM 1076 AUTPAK(I,J) = APAK(IJ+I) 1077 AUTPAK(J,I) = SIGN*APAK(IJ+I) 1078 150 CONTINUE 1079 IJ = IJ + NDIM-J 1080 200 CONTINUE 1081 END IF 1082C 1083 NTEST = 0 1084 IF( NTEST .NE. 0 ) THEN 1085 WRITE(LUPRI,*) ' AUTPAK AND APAK FROM TRIPK2 ' 1086 CALL WRTMAT(AUTPAK,NDIM,MATDIM,NDIM,MATDIM,0) 1087 CALL PRSM2(APAK,NDIM) 1088 END IF 1089C 1090 RETURN 1091 END 1092C /* Deck tribl3 */ 1093 SUBROUTINE TRIBL3(AI,AO,NRC,NCC,LRC,LCC,IBLKI,IWAY, 1094 & IBLKO,IDOBLK,SIGN) 1095C 1096C TRANSFER BETWEEN LOWER PACKED FORM AND COMPLETE FORM OF 1097C A SYMMETRIC OR ANTISYMMETRIC BLOCKED MATRIX . 1098C 1099C IWAY = 1 : UNPACKED TO PACKED FORM 1100C LOWER HALF OF COMPLETE MATRIX IS PACKED 1101C 1102C IWAY = 2 : PACKED TO UNPACKED FORM 1103C LOWER HALF OF COMPLETE MATRIX EQUAL PACKED FORM 1104C UPPER HALF IS MULTIPLIED WITH SIGN 1105C IF IDOBLK DIFFERS FROM ZERO IBLK MATRIX IS CONSTRUCTED FOR 1106C OUTPUT FORM 1107C 1108C 1109C ... JEPPE OLSEN JANUARY 1989 1110C 1111#include "implicit.h" 1112#include "priunit.h" 1113 DIMENSION AI(*),AO(*) 1114 DIMENSION LRC(NRC),LCC(NCC) 1115 DIMENSION IBLKI(NRC,NCC),IBLKO(NRC,NCC) 1116C 1117 NTEST = 0 1118C? WRITE(LUPRI,*) ' ENTERING TRIBL3 ' 1119 IF( IWAY .EQ. 1 ) THEN 1120 IF(IDOBLK.NE. 0 ) CALL ISETVC(IBLKO,0,NRC*NCC) 1121 IOFFO = 1 1122 DO 10 IRC = 1, NRC 1123 DO 5 ICC = 1, IRC 1124C? WRITE(LUPRI,*) ' IRC ICC IBKLI ', IRC,ICC,IBLKI(IRC,ICC) 1125 IF(IBLKI(IRC,ICC).GT. 0 ) THEN 1126 IF(IDOBLK .NE. 0 ) IBLKO(IRC,ICC) = IOFFO 1127 IF(IRC.EQ.ICC) THEN 1128 CALL TRIPK2(AI(IBLKI(IRC,ICC)),AO(IOFFO),1, 1129 & LRC(IRC),LCC(ICC),1.0D0) 1130 IOFFO = IOFFO + LRC(IRC)*(LRC(IRC)+1) / 2 1131 ELSE 1132 CALL COPVEC(AI(IBLKI(IRC,ICC)),AO(IOFFO), 1133 & LRC(IRC)*LCC(ICC) ) 1134 IOFFO = IOFFO + LRC(IRC)*LCC(ICC) 1135 END IF 1136 END IF 1137 5 CONTINUE 1138 10 CONTINUE 1139C 1140 IF( NTEST .NE. 0 ) THEN 1141 WRITE(LUPRI,*) 1142 & ' TRIBL3 : UNPACKED MATRIX(INPUT), PACKED MATRIX(OUTPUT)' 1143 CALL PRBL3(AI,NRC,LRC,NCC,LCC,IBLKI) 1144 WRITE(LUPRI,*) 1145C PRSBL3(A,NRC,LRC,NCC,LCC,IBLK) 1146 CALL PRSBL3(AO,NRC,LRC,NCC,LCC,IBLKO) 1147 END IF 1148 ELSE IF ( IWAY .EQ. 2 ) THEN 1149 IOFFO = 1 1150 DO 20 IRC = 1, NRC 1151 DO 15 ICC = 1, NCC 1152 IF(IRC.GE.ICC .AND. IBLKI(IRC,ICC) .GT. 0 .OR. 1153 & IRC.LT.ICC .AND. IBLKI(ICC,IRC) .GT. 0 ) THEN 1154 IF(IDOBLK.NE.0 ) IBLKO(IRC,ICC) = IOFFO 1155 IF( IRC.GT.ICC) THEN 1156 IOFFI = IBLKI(IRC,ICC) 1157 CALL COPVEC(AI(IOFFI),AO(IOFFO),LRC(IRC)*LCC(ICC)) 1158 ELSE IF ( IRC .EQ. ICC ) THEN 1159 IOFFI = IBLKI(IRC,ICC) 1160 CALL TRIPK2(AO(IOFFO),AI(IOFFI),2,LRC(IRC),LRC(IRC), 1161 & SIGN) 1162 ELSE IF ( IRC .LT. ICC) THEN 1163 IOFFI = IBLKI(ICC,IRC) 1164 CALL TRPMAT(AI(IOFFI),LCC(ICC),LRC(IRC),AO(IOFFO) ) 1165 IF(SIGN .EQ. -1.0D0) 1166 & CALL SCALVE(AO(IOFFO),SIGN,LRC(IRC)*LCC(ICC) ) 1167 END IF 1168 IOFFO = IOFFO + LRC(IRC)*LCC(ICC) 1169 END IF 1170 15 CONTINUE 1171 20 CONTINUE 1172 IF( NTEST .NE. 0 ) THEN 1173 WRITE(LUPRI,*) 1174 & ' TRIBL3 : UNPACKED MATRIX(OUTPUT), PACKED MATRIX(INPUT)' 1175 CALL PRBL3(AO,NRC,LRC,NCC,LCC,IBLKO) 1176 WRITE(LUPRI,*) 1177C PRSBL3(A,NRC,LRC,NCC,LCC,IBLK) 1178 CALL PRSBL3(AI,NRC,LRC,NCC,LCC,IBLKI) 1179 END IF 1180 END IF 1181C 1182 RETURN 1183 END 1184C /* Deck cdiag */ 1185 SUBROUTINE CDIAG(AIN,AOUT,NR,NC,FACTOR,IFLAG,ISYM) 1186C 1187C MODIFY DIAGONAL ELEMENTS OF A MATRIX 1188C 1189C IFLAG = 1 : MULTIPLY DIAGONAL ELEMENTS WITH FACTOR 1190C IFLAG = 2 : ADD DIAGONAL ELEMENTS TO FACTOR 1191C 1192C ISYM = 0 : DIAGONAL BLOCKS ARE STORED IN COMPLETE FORM 1193C ISYM .GT. 0 : DIAGONAL BLOCKS ARE STORED IN PACKED FORM 1194C ( LOWER TRIANGULAR FORM;COLUMNWISE PACKED !! ) 1195#include "implicit.h" 1196 DIMENSION AIN(*),AOUT(*) 1197C 1198 IF( ISYM .EQ. 0 ) THEN 1199 IMAX = MIN(NR,NC) 1200 CALL COPVEC(AIN,AOUT,NR*NC) 1201 IF( IFLAG .EQ. 1 ) THEN 1202 DO 100 I = 1, IMAX 1203 AOUT((I-1)*NC+I) = AIN((I-1)*NC+I) * FACTOR 1204 100 CONTINUE 1205 ELSE 1206 DO 200 I = 1, IMAX 1207 AOUT((I-1)*NC+I) = AIN((I-1)*NC+I) + FACTOR 1208 200 CONTINUE 1209 END IF 1210 ELSE IF ( ISYM .GT. 0 ) THEN 1211C NR = NC FOR SYMMETRIC PACKING 1212 NDIM = NR 1213 CALL COPVEC(AIN,AOUT,NDIM*(NDIM+1)/2) 1214 IF( IFLAG .EQ. 1 ) THEN 1215 DO 1100 I = 1, NDIM 1216 II = (I-1)*NDIM -I*(I-1)/2 + I 1217 AOUT(II) = AIN(II) * FACTOR 1218 1100 CONTINUE 1219 ELSE 1220 DO 1200 I = 1, NDIM 1221 II = (I-1)*NDIM -I*(I-1)/2 + I 1222 AOUT(II) = AIN(II) + FACTOR 1223 1200 CONTINUE 1224 END IF 1225 END IF 1226C 1227 RETURN 1228 END 1229C /* Deck cdibl3 */ 1230 SUBROUTINE CDIBL3(AI,NRC,NCC,LRC,LCC,IBLK, 1231 & AO,FACTOR,IFLAG,ISYM) 1232C 1233C CHANGE DIAGONAL ELEMENTS OF A BLOCKED MATRIX 1234C IFLAG = 1 : MULTIPLY DIAGONAL ELEMENTS WITH FACTOR 1235C IFLAG = 2 : ADD DIAGONAL ELEMENTS WITH FACTOR 1236C 1237C ISYM= 0 : DIAGONAL BLOCKS STORED IN COMPLETE FORM 1238C ISYM.GT.0 : DIAGONAL BLOCKS ARE PACKED IN LOWER TRIANGULAR FORM 1239C 1240#include "implicit.h" 1241#include "priunit.h" 1242 DIMENSION AI(*),AO(*) 1243 DIMENSION LRC(NRC),LCC(NCC) 1244 DIMENSION IBLK(NRC,NCC) 1245C 1246 DO 100 IR = 1, NRC 1247 DO 100 IC = 1, NCC 1248 IF(IBLK(IR,IC) .GT. 0 ) THEN 1249 IF( IR .NE. IC ) THEN 1250 CALL COPVEC(AI(IBLK(IR,IC)),AO(IBLK(IR,IC)), 1251 & LRC(IR)*LCC(IC) ) 1252 ELSE 1253C CDIAG(AIN,AOUT,NR,NC,FACTOR,IFLAG) 1254 CALL CDIAG(AI(IBLK(IR,IC)),AO(IBLK(IR,IC)), 1255 & LRC(IR),LCC(IC),FACTOR,IFLAG,ISYM) 1256 END IF 1257 END IF 1258 100 CONTINUE 1259C 1260 NTEST = 0 1261 IF( NTEST .NE. 0 ) THEN 1262 IF( IFLAG.EQ.1 ) WRITE(LUPRI,*) 1263 & ' DIAGONAL ELEMENTS MULTIPLIES WITH ', FACTOR 1264 IF( IFLAG.EQ.2 ) WRITE(LUPRI,*) 1265 & ' DIAGONAL ELEMENTS ADDED TO ', FACTOR 1266 WRITE(LUPRI,*) ' CDIABL3 : INPUT AND OUTPUT MATRICES ' 1267C PRBL3(A,NRC,LRC,NCC,LCC,IBLK) 1268 IF( ISYM .EQ. 0 ) THEN 1269 CALL PRBL3(AI,NRC,LRC,NCC,LCC,IBLK) 1270 WRITE(LUPRI,*) 1271 CALL PRBL3(AO,NRC,LRC,NCC,LCC,IBLK) 1272 ELSEIF( ISYM .GT. 0 ) THEN 1273 CALL PRBL3(AI,NRC,LRC,NCC,LCC,IBLK) 1274 WRITE(LUPRI,*) 1275 CALL PRBL3(AO,NRC,LRC,NCC,LCC,IBLK) 1276 END IF 1277 END IF 1278C 1279 RETURN 1280 END 1281#ifdef UNDEF 1282/* Comdeck blk_comments */ 1283C 1284C VERY OFTEN QUANTUM CHEMISTRY CODES INVOLVED A LOT 1285C OF MANIPULATION OF BLOCKED MATRICES.BLOCKING OF MATRICES 1286C CAN FOR EXAMPLE RESULT FROM SYMMETRY DIVISION OR DIVISION 1287C INTO DIFFERENT OCCUPATION TYPES.IT WOULD BE CONVENIENT TO 1288C HAVE A RATHER GENERAL SET OF ROUTINES FOR BLOCKED MATRICES. 1289C THE ROUTINES IN THIS FILE IS A BEGINNING TO A START OF A 1290C FIRST ATTEMPT ON DEVELOPING THIS SUITE OF CODES . 1291C 1292C A MATRIX IS BLOCKED BY DIVIDING THE ROWS AND COLUMNS INTO 1293C SEVERAL CLASSES.ALL ELEMENTS OF A CLASS ARE CONSECUTIVE ELEMENTS. 1294C A BLOCKING IS DEFINED BY 1295C ======================== 1296C 1297C ROWS : NRC : NUMBER OF ROW CLASSES 1298C IBRC(I): ABSOLUTE PLACE OF FIRST ELEMNT IN ROW CLASS I 1299C LRC(I) : NUMBER OF ELEMENTS IN ROW CLASS I 1300C 1301C COLUMNS : NCC : NUMBER OF COLUMN CLASSES 1302C IBCC(I): ABSOLUTE PLACE OF FIRST ELEMNT IN COLUMN CLASS I 1303C LCC(I) : NUMBER OF ELEMENTS IN COLUMN CLASS I 1304C 1305C A BLOCKED MATRIX IS DEFINED BY A MATRIX IBLK(NRC,NCC), 1306C============================ 1307C IBLK(IR,IC) GT 0 : BLOCK IR,IC IS INCLUDED. 1308C FIRST ELEMNTS OF BLOCK IS IBLK(IR,IC) 1309C ELSE : BLOCK IR,IC IS NOT INCLUDED 1310C 1311C A BLOCKED MATRIX IS STORED AS 1312C ============================== 1313C 1314C LOOP OVER ROW CLASSES 1315C LOOP OVER COLUMN CLASSES ALLOWED FOR THE ROW CLASS 1316C LOOP OVER COLUMNS OF BLOCK 1317C LOOP OVER ROWS OF BLOCK 1318C END OVER LOOP OVER ROWS OF BLOCK 1319C END OF LOOP OVER COLUMNS OF BLOCK 1320C END OVER LOOPS OF ALLOWED COLUMN BLOCKS FOR GIVEN ROW BLOCK 1321C END OF LOOP OVER ROW BLOCKS 1322C 1323C 1324C DEFINITION OF ATTRIBUTES FOR CLASSES AND ELEMENTS WILL SOON FOLLOW 1325C 1326C SUBROUTINES WRITTEN 1327C PRBL3 : PRINT BLOCKED MATRIX A 1328C PRBL3(A,NRC,LRC,NCC,LCC,IBLK) 1329C 1330#endif 1331C /* Deck prbl3 */ 1332 SUBROUTINE PRBL3(A,NRC,LRC,NCC,LCC,IBLK) 1333C 1334C PRINT BLOCKED MATRIX 1335C 1336C JANUARY 1989 , JEPPE OLSEN 1337C 1338#include "implicit.h" 1339#include "priunit.h" 1340C 1341 DIMENSION A(*) 1342 DIMENSION LRC(NRC),LCC(NCC),IBLK(NRC,NCC) 1343C 1344C WRITE(LUPRI,*) ' ENTERING PRBLM3 : ' 1345 DO 200 IRC = 1, NRC 1346 DO 100 ICC = 1, NCC 1347C WRITE(LUPRI,*) ' IRC ICC IBLK ', IRC,ICC,IBLK(IRC,IRC) 1348 IF( IBLK(IRC,ICC).GT.0 .AND. LRC(IRC)*LCC(ICC).NE.0) THEN 1349 WRITE(LUPRI,'(A,2I3)') 1350 & ' BLOCK...',IRC,ICC 1351 WRITE(LUPRI,'(A,2I3)') 1352 & ' ==================' 1353 WRITE(LUPRI,*) 1354 IPNTR = IBLK(IRC,ICC) 1355 CALL WRTMAT(A(IPNTR),LRC(IRC),LCC(ICC),LRC(IRC),LCC(ICC),0) 1356 WRITE(LUPRI,*) 1357 END IF 1358 100 CONTINUE 1359 200 CONTINUE 1360C 1361 RETURN 1362 END 1363C /* Deck tpbl3 */ 1364 SUBROUTINE TPBL3(AI,NRCI,NCCI,LRCI,LCCI,IBRCI,IBCCI,IBLKI, 1365 & AO,NRCO,NCCO,LRCO,LCCO,IBRCO,IBCCO,IBLKO ) 1366C 1367C TRANSPOSE BLOCKED MATRIX AI AND STORE IN AO 1368C ARRAYS FOR TRANSPOSED BLOCKING IS ALSE OBTAINED. 1369C 1370#include "implicit.h" 1371#include "priunit.h" 1372 DIMENSION AI(*),AO(*) 1373 DIMENSION LRCI(NRCI),LCCI(NCCI),IBRCI(NRCI),IBCCI(NCCI) 1374 DIMENSION LRCO(NCCI),LCCO(NRCI),IBRCO(NCCI),IBCCO(NRCI) 1375 DIMENSION IBLKI(NRCI,NCCI),IBLKO(NCCI,NRCI) 1376C 1377C SET UP ARRAYS FOR TRANSPOSED BLOCKING 1378C 1379 NRCO = NCCI 1380 NCCO = NRCI 1381 WRITE(LUPRI,*) ' NRCO AND NCCO ', NRCO,NCCO 1382C 1383 CALL ICOPVE(LCCI,LRCO,NRCO) 1384 CALL ICOPVE(LRCI,LCCO,NCCO) 1385 CALL ICOPVE(IBCCI,IBRCO,NRCO) 1386 CALL ICOPVE(IBRCI,IBCCO,NCCO) 1387 WRITE(LUPRI,*) ' LRCO LCCO IBRCO IBCCO ' 1388 CALL IWRTMA(LRCO,1,NRCO,1,NRCO) 1389 CALL IWRTMA(LCCO,1,NCCO,1,NCCO) 1390 CALL IWRTMA(IBRCO,1,NRCO,1,NRCO) 1391 CALL IWRTMA(IBCCO,1,NCCO,1,NCCO) 1392C 1393C TRANSPOSE FILE AND GENERATE IBLKO 1394C 1395 CALL ISETVC(IBLKO,0,NRCO*NCCO) 1396 IPNTR = 1 1397 DO 200 IRO = 1, NRCO 1398 DO 100 ICO = 1, NCCO 1399 IF(IBLKI(ICO,IRO) .GT. 0 ) THEN 1400 CALL TRPMAT(AI(IBLKI(ICO,IRO)),LRCO(ICO),LCCO(IRO), 1401 & AO(IPNTR)) 1402 IBLKO(IRO,ICO) = IPNTR 1403 IPNTR = IPNTR + LRCO(IRO)*LCCO(ICO) 1404 END IF 1405 100 CONTINUE 1406 200 CONTINUE 1407C 1408 NTEST = 1 1409 IF( NTEST .GE. 1 ) THEN 1410 WRITE(LUPRI,*) ' INFO FROM TPBL3 : ' 1411 WRITE(LUPRI,*) ' IBKLO ARRAY ' 1412 CALL IWRTMA(IBLKO,NRCO,NCCO,NRCO,NCCO) 1413 WRITE(LUPRI,*) ' INPUT MATRIX ' 1414 CALL PRBL3(AI,NRCI,LRCI,NCCI,LCCI,IBLKI) 1415 WRITE(LUPRI,*) ' OUTPUT MATRIX ' 1416 CALL PRBL3(AO,NRCO,LRCO,NCCO,LCCO,IBLKO) 1417C PRBL3(A,NRC,LRC,NCC,LRC,IBLK) 1418 END IF 1419C 1420 RETURN 1421 END 1422 FUNCTION H2TUVX(I,J,K,L,H2AC,ILTSOB) 1423C 1424C 29-Jan-1989 hjaaj 1425C l.r. 18-apr-1990 hjaaj 1426C 1427C get (tu|vx) integral where t,u,v,x are in Lunar order 1428C 1429#include "implicit.h" 1430 DIMENSION H2AC(NNASHX,NNASHX), ILTSOB(NASHT) 1431#include "iratdef.h" 1432C 1433C Used from common blocks: 1434C INFORB : NASHT,NNASHX 1435C INFTAP : LUH2AC 1436C CBGETDIS : IADH2 1437C 1438#include "inforb.h" 1439#include "inftap.h" 1440#include "cbgetdis.h" 1441C 1442C Obtain index in Sirius order 1443C 1444 I1 = ILTSOB(I) 1445 J1 = ILTSOB(J) 1446 K1 = ILTSOB(K) 1447 L1 = ILTSOB(L) 1448 IF (I1 .LT. J1) THEN 1449 ISWP = I1 1450 I1 = J1 1451 J1 = ISWP 1452 END IF 1453 IF (K1 .LT. L1) THEN 1454 ISWP = K1 1455 K1 = L1 1456 L1 = ISWP 1457 END IF 1458 IJ1 = (I1*I1 - I1)/2 + J1 1459 KL1 = (K1*K1 - K1)/2 + L1 1460C 1461C Now we can pick up the desired integral 1462C 1463 IF (IADH2 .GE. 0) THEN 1464 CALL READDX(LUH2AC,IADH2+KL1,IRAT*IJ1,H2AC) 1465 KL1 = 1 1466 END IF 1467 H2TUVX = H2AC(IJ1,KL1) 1468 RETURN 1469 END 1470C /* Deck eigen */ 1471 SUBROUTINE EIGEN(A,R,N,MV,MFKR) 1472#include "implicit.h" 1473#include "priunit.h" 1474 DIMENSION A(*),R(*) 1475 DATA TESTIT/1.D-20/ 1476 DATA TESTX /1.D-26/ 1477 DATA TESTY /1.D-18/ 1478C 1479C PURPOSE 1480C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC 1481C MATRIX 1482C 1483C USAGE 1484C CALL EIGEN(A,R,N,MV,MFKR) 1485C 1486C DESCRIPTION OF PARAMETERS 1487C A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. 1488C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF 1489C MATRIX A IN ASSCENDING ORDER. 1490C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, 1491C IN SAME SEQUENCE AS EIGENVALUES) 1492C N - ORDER OF MATRICES A AND R 1493C MV- INPUT CODE 1494C 0 COMPUTE EIGENVALUES AND EIGENVECTORS 1495C 1 COMPUTE EIGENVALUES ONLY (R NEED NOT BE 1496C DIMENSIONED BUT MUST STILL APPEAR IN CALLING 1497C SEQUENCE) 1498C MFKR=0 NO SORT 1499C =1 SORT 1500C 1501C REMARKS 1502C ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1) 1503C MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R 1504C 1505C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 1506C NONE 1507C 1508C METHOD 1509C DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED 1510C BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN "MATHEMATICAL 1511C METHODS FOR DIGITAL COMPUTERS", EDITED BY A. RALSTON AND 1512C H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7 1513C 1514C .................................................................. 1515C 1516C 1517C ............................................................... 1518C 1519C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE 1520C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION 1521C STATEMENT WHICH FOLLOWS. 1522C 1523C DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX, 1524C 1 COSX2,SINCS,RANGE 1525C 1526C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS 1527C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS 1528C ROUTINE. 1529C 1530C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO 1531C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. SQRT IN STATEMENTS 1532C 40, 68, 75, AND 78 MUST BE CHANGED TO DSQRT. ABS IN STATEMENT 1533C 62 MUST BE CHANGED TO DABS. THE CONSTANT IN STATEMENT 5 SHOULD 1534C BE CHANGED TO 1.0D-12. 1535C 900111-hjaaj: use generic SQRT and ABS 1536C 1537C ............................................................... 1538C 1539C GENERATE IDENTITY MATRIX 1540C 1541 RANGE=1.0D-12 1542 IF(MV.eq.1) GO TO 25 1543 IQ=-N 1544 DO 20 J=1,N 1545 IQ=IQ+N 1546 DO 20 I=1,N 1547 IJ=IQ+I 1548 R(IJ)=0.0D+00 1549 IF(I.eq.J) R(IJ)=1.0D+00 1550 20 CONTINUE 1551C 1552C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) 1553C 1554 25 ANORM=0.0D+00 1555 DO 35 J=2,N 1556 DO 35 I=1,J-1 1557 IA=I+(J*J-J)/2 1558 ANORM=ANORM+A(IA)*A(IA) 1559 35 CONTINUE 1560 IF(ANORM .LE. 0.0D0) GO TO 165 1561 ANORM=SQRT(2.0D0*ANORM) 1562 ANRMX=ANORM*RANGE/DFLOAT(N) 1563C 1564C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR 1565C 1566 IND=0 1567 THR=ANORM 1568 45 THR=THR/DFLOAT(N) 1569 IF(THR.LT.TESTY)THR=0.D0 1570 50 L=1 1571 55 M=L+1 1572C 1573C COMPUTE SIN AND COS 1574C 1575 60 MQ=(M*M-M)/2 1576 LQ=(L*L-L)/2 1577 LM=L+MQ 1578 IF(ABS(A(LM)).LT.TESTY)A(LM)=0.D0 1579 IF(ABS(A(LM)).EQ.0.D0.AND.THR.EQ.0.D0)GO TO 130 1580 IF(ABS(A(LM)) .lt. THR) go to 130 1581 IND=1 1582 LL=L+LQ 1583 MM=M+MQ 1584 X=0.5D+00*(A(LL)-A(MM)) 1585 AJUK=(A(LM)*A(LM)+X*X) 1586 AJUK=SQRT(AJUK) 1587 IF(ABS(AJUK).LT.TESTIT) THEN 1588 WRITE(LUPRI,3000)TESTIT,AJUK,A(LM) 1589 3000 FORMAT(' ***in EIGEN: DENOMINATOR LT ',D14.6,'. VALUE=',D14.6, 1590 &'. NUMERATOR=',D14.6) 1591 Y=0.D0 1592 ELSE 1593 Y=-A(LM)/AJUK 1594 END IF 1595C Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X) 1596 IF(X .lt. 0.0d0) Y=-Y 1597 AJUK=(1.D0-Y*Y) 1598 IF(AJUK.LT.0.D0) THEN 1599 WRITE(LUPRI,3001) AJUK 1600 3001 FORMAT(' ***in EIGEN: SQRT OF ',D20.8) 1601 AJUK=0.D0 1602 ELSE 1603 AJUK=SQRT(AJUK) 1604 END IF 1605 AJUK=2.D0*(1.D0+AJUK) 1606 AJUK=SQRT(AJUK) 1607 SINX=Y/AJUK 1608C SINX=Y/ SQRT(2.0D+00*(1.0D+00+( SQRT(1.0D+00-Y*Y)))) 1609 SINX2=SINX*SINX 1610C COSX= SQRT(1.0D+00-SINX2) 1611 AJUK=1.D0-SINX2 1612 IF(AJUK.LT.TESTX)AJUK=0.D0 1613 COSX=SQRT(AJUK) 1614 COSX2=COSX*COSX 1615 SINCS =SINX*COSX 1616C 1617C ROTATE L AND M COLUMNS 1618C 1619 ILQ=N*(L-1) 1620 IMQ=N*(M-1) 1621 DO 125 I=1,N 1622 IQ=(I*I-I)/2 1623 IF(I.eq.L .or. I.eq.M) go to 115 1624 IF(I.lt.M) THEN 1625 IM=I+MQ 1626 ELSE 1627 IM=M+IQ 1628 END IF 1629 IF(I.lt.L) THEN 1630 IL=I+LQ 1631 ELSE 1632 IL=L+IQ 1633 END IF 1634 110 X=A(IL)*COSX-A(IM)*SINX 1635 A(IM)=A(IL)*SINX+A(IM)*COSX 1636 A(IL)=X 1637 115 IF(MV.eq.1) go to 125 1638 ILR=ILQ+I 1639 IMR=IMQ+I 1640 X=R(ILR)*COSX-R(IMR)*SINX 1641 R(IMR)=R(ILR)*SINX+R(IMR)*COSX 1642 R(ILR)=X 1643 125 CONTINUE 1644 X=2.0D+00*A(LM)*SINCS 1645 Y=A(LL)*COSX2+A(MM)*SINX2-X 1646 X=A(LL)*SINX2+A(MM)*COSX2+X 1647 A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) 1648 A(LL)=Y 1649 A(MM)=X 1650C 1651C TESTS FOR COMPLETION 1652C 1653C TEST FOR M = LAST COLUMN 1654C 1655 130 IF (M.ne.N) THEN 1656 M=M+1 1657 GO TO 60 1658 END IF 1659C 1660C TEST FOR L = SECOND FROM LAST COLUMN 1661C 1662 IF (L.ne.(N-1)) THEN 1663 L=L+1 1664 GO TO 55 1665 END IF 1666 IF (IND.eq.1) THEN 1667 IND=0 1668 GO TO 50 1669 END IF 1670C 1671C COMPARE THRESHOLD WITH FINAL NORM 1672C 1673 IF(THR .gt. ANRMX) go to 45 1674C 1675C SORT EIGENVALUES AND EIGENVECTORS 1676C 1677 165 IQ=-N 1678 IF(MFKR.EQ.0)GO TO 186 1679 DO 185 I=1,N 1680 IQ=IQ+N 1681 LL=I+(I*I-I)/2 1682 JQ=N*(I-2) 1683 DO 185 J=I,N 1684 JQ=JQ+N 1685 MM=J+(J*J-J)/2 1686 IF(A(MM).ge.A(LL)) go to 185 1687 170 X=A(LL) 1688 A(LL)=A(MM) 1689 A(MM)=X 1690 IF(MV.eq.1) go to 185 1691 DO 180 K=1,N 1692 ILR=IQ+K 1693 IMR=JQ+K 1694 X=R(ILR) 1695 R(ILR)=R(IMR) 1696 180 R(IMR)=X 1697 185 CONTINUE 1698186 CONTINUE 1699 RETURN 1700 END 1701