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 flshfo */ 20 SUBROUTINE FLSHFO (IUNIT) 21C 22C *** THIS SUBROUTINE IS SYSTEM DEPENDENT *** 23C 24C Flush formatted output unit (empty buffers). 25C If no flush utility, this is achieved by 26C CLose and reOPen Formatted Output 27C 28C Written 21-Nov-1983 by Hans Jorgen Aa. Jensen in Uppsala, Sweden. 29C Last revision 16-Jul-1984 hjaaj / 30-Oct-1984 hjaaj (extendsize) 30C 10-Feb-1989 hjaaj, renamed CLOPFO to FLSHFO 31C 32C Calls to this subroutine makes it possible to read the output 33C up to the moment of the last call while the program continues 34C executing (provided the computer allows shared access). 35C This subroutine may be a dummy routine. 36C 37 IF (IUNIT .GE. 0) THEN 38C ... do not try to flush unused units (e.g. LUPRI on a slave) /hjaaj 39#if defined (SYS_AIX) || defined (SYS_LINUX) || defined (SYS_DARWIN) 40C 41C Force transfer of all buffered output to the file or device 42C associated with logical unit IUNIT. 43C 44 CALL FLUSH(IUNIT) 45#endif 46 END IF 47 RETURN 48 END 49C /* Deck getdat */ 50 SUBROUTINE GETDAT(CDATE,CTIME) 51C 52C 24-Jan-1988 Hans Joergen Aa. Jensen 53C 54C Return date and time as character*8, for labels. 55C 56 CHARACTER*(8) CDATE, CTIME 57C ... do not try to flush unused units (e.g. LUPRI on a slave) /hjaaj 58#if defined (SYS_AIX) || defined (SYS_LINUX) || defined (SYS_DARWIN) 59 CHARACTER*(24) FFDATE 60 CHARACTER*(24) FDATE 61 FFDATE = FDATE() 62 CDATE = FFDATE(9:10)//FFDATE(5:7)//FFDATE(23:24)//' ' 63 CTIME = FFDATE(12:19) 64#else 65 CDATE = ' -date- ' 66 CTIME = ' -time- ' 67#endif 68 RETURN 69 END 70 71 subroutine gettim(cputime,walltime) 72 ! returns used CPU time and elapsed wall time 73 ! since first call of GETTIM 74 implicit none 75 real*8 cputime, walltime 76 77 logical first 78 data first /.true./ 79 80 real*8 TCPU0, TWALL0, tcpu1, twall1 81 save TCPU0, TWALL0 82 integer*4 dateandtime0(8), dateandtime1(8) 83 84 if (first) then 85 first = .false. 86 call cpu_time(TCPU0) 87 call date_and_time(values=dateandtime0) 88 call get_walltime(dateandtime0,TWALL0) 89 end if 90 call cpu_time(tcpu1) 91 call date_and_time(values=dateandtime1) 92 call get_walltime(dateandtime1,twall1) 93 94 cputime = tcpu1 - TCPU0 95 walltime = twall1 - TWALL0 96 97 end subroutine gettim 98 99!> \brief Get elapsed walltime in seconds since 1/1-2010 00:00:00 100!> \author S. Host 101!> \date October 2010 102!> 103!> Years that are evenly divisible by 4 are leap years. 104!> Exception: Years that are evenly divisible by 100 are not leap years, 105!> unless they are also evenly divisible by 400. Source: Wikipedia 106!> 107 subroutine get_walltime(dateandtime,walltime) 108 implicit none 109 integer*4 dateandtime(8) ! "values" output from fortran intrinsic subroutine date_and_time 110 real*8 walltime ! Elapsed wall time in seconds 111 integer*4 month, year 112 113! The output from the fortran intrinsic date_and_time 114! gives the following values: 115! 1. Year 116! 2. Month 117! 3. Day of the month 118! 4. Time difference in minutes from Greenwich Mean Time (GMT) 119! 5. Hour 120! 6. Minute 121! 7. Second 122! 8. Millisecond 123 124! Count seconds, minutes, hours, days, months and years and sum up seconds: 125 126 walltime = 1.0d-3*dateandtime(8) !Milliseconds 127 walltime = walltime + dateandtime(7)*1.0d0 !Seconds counted 128 walltime = walltime + 60d0*dateandtime(6) !Minutes counted 129 walltime = walltime + 3600d0*dateandtime(5) !Hours counted 130 walltime = walltime + 24d0*3600d0*(dateandtime(3)-1) !Days counted (substract 1 to count only whole days) 131 132 !Months are special, since they are not equally long: 133 134 do month = 1, dateandtime(2)-1 !substract 1 to count only whole months 135 if (month == 1 .or. month == 3 .or. month == 5 .or. 136 & month == 7 .or. month == 8 .or. month == 10) then !Since we subtract 1, month can never be 12 137 walltime = walltime + 31d0*24d0*3600d0 138 else if (month == 2) then 139 if (.false.) then !insert exception for if current year is a leap year 140 walltime = walltime + 29d0*24d0*3600d0 141 else 142 walltime = walltime + 28d0*24d0*3600d0 143 endif 144 else if (month == 4 .or. month == 6 .or. month == 9 .or. 145 & month == 11) then 146 walltime = walltime + 30d0*24d0*3600d0 147 else 148 call quit('Unknown month (get_walltime)') 149 endif 150 enddo 151 152 !Years are special, since leap years are one day longer: 153 154 do year = 2010, dateandtime(1) 155 if (mod(year,400)==0) then 156 walltime = walltime + 366*24*3600 !Leap year 157 else if (mod(year,100)==0) then 158 walltime = walltime + 365*24*3600 !Not leap year 159 else if (mod(year,4)==0) then 160 walltime = walltime + 366*24*3600 !Leap year 161 else 162 walltime = walltime + 365*24*3600 !Not leap year 163 endif 164 enddo 165 166 end subroutine get_walltime 167C /* Deck timtxt */ 168 SUBROUTINE TIMTXT(TEXT,TIMUSD,LUPRIN) 169C 170C TIMTXT based on TIMER by TUH //900709-hjaaj 171C 172#include "implicit.h" 173 CHARACTER*(*) TEXT 174 CHARACTER AHOUR*6, ASEC*8, AMIN*8 175C 176 ISECND = NINT(TIMUSD) 177 IF (ISECND .GE. 60) THEN 178 MINUTE = ISECND/60 179 IHOURS = MINUTE/60 180 MINUTE = MINUTE - 60*IHOURS 181 ISECND = ISECND - 3600*IHOURS - 60*MINUTE 182 IF (IHOURS .EQ. 1) THEN 183 AHOUR = ' hour ' 184 ELSE 185 AHOUR = ' hours' 186 END IF 187 IF (MINUTE .EQ. 1) THEN 188 AMIN = ' minute ' 189 ELSE 190 AMIN = ' minutes' 191 END IF 192 IF (ISECND .EQ. 1) THEN 193 ASEC = ' second ' 194 ELSE 195 ASEC = ' seconds' 196 END IF 197 IF (IHOURS .GT. 0) THEN 198 WRITE(LUPRIN,100) 199 * TEXT, IHOURS, AHOUR, MINUTE, AMIN, ISECND, ASEC 200 ELSE 201 WRITE(LUPRIN,200) TEXT, MINUTE, AMIN, ISECND, ASEC 202 END IF 203 ELSE 204 WRITE(LUPRIN,300) TEXT,TIMUSD 205 END IF 206 100 FORMAT(1X,A,I4,A,I3,A,I3,A) 207 200 FORMAT(1X,A, I3,A,I3,A) 208 300 FORMAT(1X,A,F7.2,' seconds') 209 RETURN 210 END 211C /* Deck tstamp */ 212 SUBROUTINE TSTAMP(TEXT,LUPRIN) 213C 214C Copyright Hans Joergen Aa. Jensen 9-Jul-1990 215C 216C Purpose: To stamp as many as possible of 217C text, date, time, computer, and hostname to LUPRIN 218C 219#include "implicit.h" 220 CHARACTER*(*) TEXT 221C 222#if defined (SYS_UNIX) || defined (SYS_DARWIN) || defined (SYS_LINUX) 223 CHARACTER*(40) HSTNAM 224 CHARACTER*(24) FDATE 225#endif 226#if defined (SYS_AIX) 227 CHARACTER*(24) fdate 228 CHARACTER*(32) HSTNAM 229#endif 230C 231 LTEXT = LEN(TEXT) 232 IF (LTEXT .GT. 0) THEN 233 IF (TEXT(1:LTEXT) .EQ. 'INIT') THEN 234CHJ March 2005: this is not working on all machines, so ... 235C WRITE (LUPRIN,'(3(/T3,2A)/)') 236C & 'Last compilation : ', 237C & LAST_DALTON_COMPILATION 238C & ,'Fortran and C compilers: ', 239C & F_AND_C_COMPILERS 240C & ,'Scientific libraries : ', 241C & LIBRARY_LIST 242 WRITE (LUPRIN,'()') 243 ELSE 244 WRITE (LUPRIN,'(/A)') TEXT 245 END IF 246 ELSE 247 WRITE (LUPRIN,'()') 248 END IF 249#if defined (SYS_LINUX) 250 WRITE (LUPRIN,'(T6,2A)') 'Date and time (Linux) : ',FDATE() 251#endif 252#if defined (SYS_DARWIN) 253 WRITE (LUPRIN,'(T6,2A)') 'Date and time (Darwin) : ',FDATE() 254#endif 255#if defined (SYS_AIX) 256 WRITE (LUPRIN,'(T6,2A)') 'Date and time (IBM-AIX): ',fdate() 257#endif 258 259#if defined (SYS_LINUX) || defined (SYS_UNIX) || defined (SYS_AIX) || defined (SYS_DARWIN) 260 CALL HOSTNM(HSTNAM) 261 WRITE (LUPRIN,'(T6,2A)') 'Host name : ',HSTNAM 262#endif 263C 264 RETURN 265 END 266C /* Deck ordrss */ 267 SUBROUTINE ORDRSS(EVEC,EVAL,ISS,N,NEVEC) 268C 269C 920729-hjaaj (based on ORDER) 270C 271C Purpose: order the N values in EVAL and their associated vectors 272C in EVEC so EVAL(i+1) .ge. EVAL(i), 273C but only within the class of vectors having the 274C same value in the ISS array (which could be the 275C supersymmetry of the orbital). 276C 277#include "implicit.h" 278 DIMENSION EVEC(*),EVAL(*),ISS(*) 279 IF (N.LE.1) RETURN 280 IN = 1 281 DO 10 I=1,N-1 282 EMIN = EVAL(I) 283 IMIN = I 284 ISSI = ISS(I) 285 DO 20 J=I+1,N 286C IF (ISS(J) .NE. ISSI) GO TO 20 287C ... now also reorder diff. supsym, instead update ISS(:) 288C /hjaaj aug 04 289 IF (EVAL(J) .LT. EMIN) THEN 290 EMIN = EVAL(J) 291 IMIN = J 292 ENDIF 293 20 CONTINUE 294 IF (IMIN.NE.I) THEN 295 EVAL(IMIN)=EVAL(I) 296 EVAL(I)=EMIN 297 IF (NEVEC .GT. 0) THEN 298 CALL DSWAP(NEVEC,EVEC(IN),1,EVEC((IMIN-1)*NEVEC+1),1) 299 ENDIF 300 ISS(I) = ISS(IMIN) 301 ISS(IMIN) = ISSI 302 ENDIF 303 IN = IN + NEVEC 304 10 CONTINUE 305 RETURN 306 END 307C /* Deck ord2ss */ 308 SUBROUTINE ORD2SS(EVEC,EVAL,ISS,N,NEVEC) 309C 310C Purpose: order the N values in EVAL and their associated vectors 311C in EVEC so EVAL(i+1) .le. EVAL(i) using the infomation in ISS 312C (this is opposite order of "ORDRSS") 313C 314#include "implicit.h" 315 DIMENSION EVEC(*),EVAL(*),ISS(*) 316 IF (N.LE.1) RETURN 317 IN = 1 318 DO 10 I=1,N-1 319 EMAX = EVAL(I) 320 IMAX = I 321 ISSI = ISS(I) 322 DO 20 J=I+1,N 323C IF (ISS(J) .NE. ISSI) GO TO 20 324C ... now also reorder diff. supsym, instead update ISS(:) 325C /hjaaj aug 04 326 IF (EVAL(J) .GT. EMAX) THEN 327 EMAX = EVAL(J) 328 IMAX = J 329 ENDIF 330 20 CONTINUE 331 IF (IMAX.NE.I) THEN 332 EVAL(IMAX)=EVAL(I) 333 EVAL(I)=EMAX 334 IF (NEVEC .GT. 0) THEN 335 CALL DSWAP(NEVEC,EVEC(IN),1,EVEC((IMAX-1)*NEVEC+1),1) 336 ENDIF 337 ISS(I) = ISS(IMAX) 338 ISS(IMAX) = ISSI 339 ENDIF 340 IN = IN + NEVEC 341 10 CONTINUE 342 RETURN 343 END 344C /* Deck order */ 345 SUBROUTINE ORDER(EVEC,EVAL,N,NEVEC) 346C 347C Purpose: order the N values in EVAL and their associated vectors 348C in EVEC so EVAL(i+1) .ge. EVAL(i) 349C 350C Revisions: 351C 29-Jul-1992 hjaaj (only dswap if nevec .gt. 0) 352C 2-Nov-1984 hjaaj (new parameter NEVEC, EVEC(1:NEVEC,1:N)) 353C 27-Oct-1984 hjaaj (reduced number of swaps) 354C 355#include "implicit.h" 356 DIMENSION EVEC(*),EVAL(*) 357 IF (N.LE.1) RETURN 358 IN = 1 359 DO 10 I=1,N-1 360 EMIN = EVAL(I) 361 IMIN = I 362 DO 20 J=I+1,N 363 IF (EVAL(J) .LT. EMIN) THEN 364 EMIN = EVAL(J) 365 IMIN = J 366 ENDIF 367 20 CONTINUE 368 IF (IMIN.NE.I) THEN 369 EVAL(IMIN)=EVAL(I) 370 EVAL(I)=EMIN 371 IF (NEVEC .GT. 0) THEN 372 CALL DSWAP(NEVEC,EVEC(IN),1,EVEC((IMIN-1)*NEVEC+1),1) 373 ENDIF 374 ENDIF 375 IN = IN + NEVEC 376 10 CONTINUE 377 RETURN 378 END 379C /* Deck order2 */ 380 SUBROUTINE ORDER2(EVEC,EVAL,N,NEVEC) 381C 382C Purpose: order the N values in EVAL and their associated vectors 383C in EVEC so EVAL(i+1) .le. EVAL(i) 384C (this is opposite order of "ORDER") 385C 386C Revisions: 387C 29-Jul-1992 hjaaj (only dswap if nevec .gt. 0) 388C 5-Aug-1985 hjaaj (first version, based on ORDER) 389C 390#include "implicit.h" 391 DIMENSION EVEC(*),EVAL(*) 392 IF (N.LE.1) RETURN 393 IN = 1 394 DO 10 I=1,N-1 395 EMAX = EVAL(I) 396 IMAX = I 397 DO 20 J=I+1,N 398 IF (EVAL(J) .GT. EMAX) THEN 399 EMAX = EVAL(J) 400 IMAX = J 401 ENDIF 402 20 CONTINUE 403 IF (IMAX.NE.I) THEN 404 EVAL(IMAX)=EVAL(I) 405 EVAL(I)=EMAX 406 IF (NEVEC .GT. 0) THEN 407 CALL DSWAP(NEVEC,EVEC(IN),1,EVEC((IMAX-1)*NEVEC+1),1) 408 ENDIF 409 ENDIF 410 IN = IN + NEVEC 411 10 CONTINUE 412 RETURN 413 END 414C /* Deck our_own_traceback */ 415 SUBROUTINE OUR_OWN_TRACEBACK 416C 417C Written 4-Dec-1983 hjaaj; last revision Jan. 2011 418C 419#include "implicit.h" 420#include "priunit.h" 421 422#if defined (SYS_AIX) 423 include 'fexcp.h' 424#endif 425 426 INTEGER A,B,C 427 SAVE A,B,C 428 DATA A/1/,B/0/,C/0/ 429 430C 920522-hjaaj -- ad hoc routine for creating traceback on IBM-AIX 431C and some other operating systems 432C Note: integer divide by zero is the only error which 433C always will cause an exception 434C 435 436#if defined (SYS_AIX) 437 call SIGNAL(SIGTRAP,xl__trce) 438#else 439 440 C = C + A/B 441 442C The print statements below will only be printed if the integer divide by zero 443C code above does not cause the program to exit. 444#endif 445 446Chjaaj apr06: 447C Use stderr on unix/linux systems, if LUPRI not defined yet. 448 IF (LUPRI .LT. 0) LUPRI = 0 449#if defined (SYS_LINUX) 450 WRITE(LUPRI,*) 'Linux has no obvious system traceback facility.' 451#endif 452#if defined (SYS_DARWIN) 453 WRITE(LUPRI,*) 'Darwin has no obvious system traceback facility.' 454#endif 455 456 RETURN 457 END 458C /* Deck canon */ 459 SUBROUTINE CANON(I,J,K,L) 460C reorder I, J, K, L to canonical 2-electron integral order 461#include "implicit.h" 462 IP=MAX(I,J) 463 JP=I+J-IP 464 KP=MAX(K,L) 465 LP=K+L-KP 466 IF (IP.GT.KP) THEN 467 I=IP 468 J=JP 469 K=KP 470 L=LP 471 ELSE 472 I=KP 473 J=LP 474 K=IP 475 L=JP 476 IF(I.NE.K)RETURN 477 IF(J.GT.L)RETURN 478 J=JP 479 L=LP 480 END IF 481 RETURN 482 END 483C /* Deck jaco */ 484 SUBROUTINE JACO (F,V,NB,NMAX,NROWV,BIG,JBIG) 485C 486C F is symmetric packed matrix of dimension NB. 487C The first block of size NMAX will be diagonalized. 488C V is for eigenvectors, only V(NROWV,NMAX) will be referenced. 489C On entry it must correspond the basis vectors corresponding to the 490C F matrix on entry, e.g. unit matrix or AO coefficients for each MO. 491C 492C Revisions: 493C 2-Nov-1984 hjaaj (new parameter NROWV such that 494C dim(V) = (NROWV,NMAX). This makes 495C it possible to solve eigenproblem 496C in a reduced basis but get the 497C eigenvectors in the original full 498C basis, e.g. less mo's than ao's) 499C 23-Feb-1989 hjaaj Note that if NROWV = 0 then only 500C eigenvalues will be calculated, 501C V matrix will not be referenced. 502C 27-Jul-1990 hjaaj Changed -CX,+SX transformation to +CX,-SX 503C transformation; probably the -CX,+SX 504C transformation was responsible for that 505C the eigenvectors easily changed sign. 506C Changed initial test on NB. Changed SD. 507C Optimized IR loop. 508C Jun-1992 ov Parameters for 0.5, 1.5, ... (for Cray) 509C 20-Jul-1992 hjaaj Changed C1,C2 to THRZER 510C 30-oct-1992 hjaaj zero f(ab) to avoid round-off errors 511C absolute conv.threshold SD=C1 512C 18-aug-2005 wmk Changed C1 to 1.D-15 513#include "implicit.h" 514#include "priunit.h" 515 DIMENSION F(*),V(NROWV,*) 516 DIMENSION BIG(*) ,JBIG(*) 517 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, ROOT2 = 0.707106781186548D0) 518 PARAMETER(DP5 = 0.5D0, D1P5 = 1.5D0, D1P375 = 1.375D0, 519 * D3P875 = 3.875D0, DP25 = 0.25D0) 520#include "thrzer.h" 521 DATA C1,C2,C3,C4,C5,C6/1.D-15,THRZER,1.D-20,1.D-14,1.D-9,1.D-5/ 522Cwas: DATA C1,C2,C3,C4,C5,C6/THRZER,THRZER,1.D-20,1.D-14,1.D-9,1.D-5/ 523Cwas: DATA C1,C2,C3,C4,C5,C6/1.D-12,1.D-12,1.D-20,1.D-14,1.D-9,1.D-5/ 524 IF (NB.LE.1 .OR. NMAX.LE.0) RETURN 525 CALL QENTER('JACO') 526 CALL GETTIM(TSTRT, WSTRT) 527 DO 190 I=1,NB 528 JBIGI=0 529 J=MIN(I-1,NMAX) 530 IF (J .GT. 0) THEN 531 II = (I*I-I)/2 532 ABIGI=D0 533 DO 18 K=1,J 534 IF (ABIGI .GE. ABS(F(II+K))) GO TO 18 535 ABIGI=ABS(F(II+K)) 536 JBIGI=K 537 18 CONTINUE 538 END IF 539 IF (JBIGI .GT. 0) THEN 540 JBIG(I) = JBIGI 541 BIG(I) = F(II+JBIGI) 542 ELSE 543 JBIG(I) = 0 544 BIG(I) = D0 545 END IF 546 190 CONTINUE 547C 548#if defined (VAR_OLDCODE) 549C 900727-hjaaj: 550C SD calculation was done in every Jacobi iteration. 551C Now the largest absolute element in F is found once and 552C the SD based on that value is used in every iteration. 553 410 SD=1.05D 00 554 DO 220 J=1,NMAX 555 DAB=ABS(F(J*(J+1)/2)) 556CHJ-861103: commented out next line, it seems to make the loop 557C meaningless (setting SD equal to J=NMAX value always!) 558C IF (SD .GT. DAB) SD=DAB 559 220 SD=MAX(SD,DAB) 560 SD=MAX(C1,C2*SD) 561#else 562C 921030-hjaaj: SD = C1 now 563 NNB = (NB*NB+NB)/2 564C SD = 1.05D0 565C DO 220 J = 1,NNB 566C SD = MAX(SD, ABS(F(J)) ) 567C 220 CONTINUE 568C SD=MAX(C1,C2*SD) 569 SD=C1 570C 571 MXITJA = 50*NNB 572 ITJACO = 0 573 410 ITJACO = ITJACO + 1 574 IF (ITJACO .GT. MXITJA) THEN 575 CALL QUIT('ERROR: JACO did not converge ...') 576 END IF 577#endif 578 T = D0 579 DO 230 I=2,NB 580 IF (T .GE. ABS(BIG(I))) GO TO 230 581 T = ABS(BIG(I)) 582 IB= I 583 230 CONTINUE 584 IF(T.LT.SD) GO TO 420 585 IA =JBIG(IB) 586 IAA=IA*(IA-1)/2 587 IBB=IB*(IB-1)/2 588 DIF=F(IAA+IA)-F(IBB+IB) 589 IF( ABS(DIF) .GT. C3) GO TO 271 590 SX=ROOT2 591 CX=ROOT2 592 GO TO 270 593 271 T2X2 =BIG(IB)/DIF 594 T2X25=T2X2*T2X2 595 IF(T2X25 .GT. C4) GO TO 240 596 CX=D1 597 SX=T2X2 598 GO TO 270 599 240 IF(T2X25 .GT. C5) GO TO 250 600 SX=T2X2*(D1 - D1P5*T2X25) 601 CX=D1 - DP5*T2X25 602 GO TO 270 603 250 IF(T2X25 . GT . C6) GO TO 260 604 CX=D1+T2X25*(T2X25*D1P375 - DP5 ) 605 SX= T2X2*(D1 + T2X25*(T2X25*D3P875 - D1P5)) 606 GO TO 270 607 260 T = DP25 / SQRT(DP25 + T2X25) 608 CX= SQRT(DP5 + T) 609 SX= SIGN( SQRT(DP5 - T),T2X2) 610 270 CONTINUE 611 DO 275 IR=1,IA 612 T = F(IAA+IR)*SX 613 F(IAA+IR)= F(IAA+IR)*CX+F(IBB+IR)*SX 614 F(IBB+IR)=-T +F(IBB+IR)*CX 615 275 CONTINUE 616 IEAA=IAA+IA 617 IEAB=IBB+IA 618 TT =F(IEAB) 619 F(IEAB)=BIG(IB) 620 IF (JBIG(IA) .EQ. 0) THEN 621 IRST = IA + 1 622 IEAR = IEAA + IA 623 IEBR = IEAB + 1 624 ELSE 625 IRST = IA 626 IEAR = IEAA 627 IEBR = IEAB 628 END IF 629 DO 390 IR = IRST,NB 630#if !defined (VAR_OLDCODE) 631 IF (IR .EQ. IA) GO TO 360 632C ... we have checked above that JBIG(IA) .ne. 0 633#else 634 IF (IR .EQ. IA) THEN 635 GO TO 360 636C ... we have checked above that JBIG(IA) .ne. 0 637C IF(JBIG(IR)) 360,380,360 638 END IF 639#endif 640 T = F(IEAR)*SX 641 F(IEAR)= F(IEAR)*CX+F(IEBR)*SX 642 F(IEBR)=-T +F(IEBR)*CX 643 T =F(IEAR) 644 IT =IA 645 IF(IR-IB) 340,310,320 646 310 F(IEAA)=F(IEAA)*CX+F(IEAB)*SX 647C 921030+hjaaj: zero f(ab) to avoid round-off errors 648C F(IEAB)= TT*CX+F(IEBR)*SX 649 F(IEAB)= D0 650 F(IEBR)= -TT*SX+F(IEBR)*CX 651 GO TO 360 652 320 IF(ABS(T) .GE. ABS(F(IEBR))) GO TO 340 653 T =F(IEBR) 654 IT =IB 655 340 IF(ABS(T) .LT. ABS(BIG(IR))) GO TO 350 656 BIG(IR) = T 657 JBIG(IR) = IT 658 GO TO 380 659 350 IF(IA .NE. JBIG(IR) .AND. IB .NE. JBIG(IR)) GO TO 380 660 360 K= IEAR - IA 661 JBIGI = 0 662 IR1=MIN (IR-1,NMAX) 663 IF (IR1 .GT. 0) THEN 664 ABIGI = D0 665 DO 370 I=1,IR1 666 IF(ABIGI .GE. ABS(F(K+I))) GO TO 370 667 ABIGI = ABS(F(K+I)) 668 JBIGI =I 669 370 CONTINUE 670 END IF 671 IF (JBIGI .GT. 0) THEN 672 JBIG(IR) = JBIGI 673 BIG(IR) = F(K+JBIGI) 674 ELSE 675 JBIG(IR) = 0 676 BIG(IR) = D0 677 END IF 678 380 CONTINUE 679 IEAR = IEAR + IR 680 IF (IR .GE. IB) THEN 681 IEBR = IEBR + IR 682 ELSE 683 IEBR = IEBR + 1 684 END IF 685 390 CONTINUE 686C 687 DO I=1,NROWV 688 VIIA = V(I,IA) 689 VIIB = V(I,IB) 690 V(I,IA) = VIIA*CX + VIIB*SX 691 V(I,IB) = -VIIA*SX + VIIB*CX 692 END DO 693 GO TO 410 694 420 CONTINUE 695c CALL GETTIM(TEND, WEND) 696c WRITE(LUPRI,'(/A,4I10/A,2F20.2)') 697c & 'JACO -- ITJACO, NB,NMAX,NROWV :',ITJACO, NB,NMAX,NROWV, 698c & 'JACO -- CPUTIME, WALLTIME :',TEND-TSTRT, WEND-WSTRT 699 CALL QEXIT('JACO') 700 RETURN 701 END 702C /* Deck jaco_thr */ 703 SUBROUTINE JACO_THR (F,V,NB,NMAX,NROWV,THR_conv) 704C 705C F is symmetric packed matrix of dimension NB. 706C The first block of size NMAX will be diagonalized. 707C V is for eigenvectors, only V(NROWV,NMAX) will be referenced. 708C On entry it must correspond the basis vectors corresponding to the 709C F matrix on entry, e.g. unit matrix or AO coefficients for each MO. 710C 711C Revisions: 712C 2-Nov-1984 hjaaj (new parameter NROWV such that 713C dim(V) = (NROWV,NMAX). This makes 714C it possible to solve eigenproblem 715C in a reduced basis but get the 716C eigenvectors in the original full 717C basis, e.g. less mo's than ao's) 718C 23-Feb-1989 hjaaj Note that if NROWV = 0 then only 719C eigenvalues will be calculated, 720C V matrix will not be referenced. 721C 27-Jul-1990 hjaaj Changed -CX,+SX transformation to +CX,-SX 722C transformation; probably the -CX,+SX 723C transformation was responsible for that 724C the eigenvectors easily changed sign. 725C Changed initial test on NB. Changed SD. 726C Optimized IR loop. 727C Jun-1992 ov Parameters for 0.5, 1.5, ... (for Cray) 728C 20-Jul-1992 hjaaj Changed C1,C2 to THRZER 729C 30-oct-1992 hjaaj zero f(ab) to avoid round-off errors 730C absolute conv.threshold SD=C1 731C 18-aug-2005 wmk Changed C1 to 1.D-15 732#include "implicit.h" 733#include "priunit.h" 734 DIMENSION F(*),V(NROWV,*) 735 736! local variables and arrays: 737 738 DIMENSION BIG(NB), JBIG(NB) ! automatic arrays 739 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, ROOT2 = 0.707106781186548D0) 740 PARAMETER(DP5 = 0.5D0, D1P5 = 1.5D0, D1P375 = 1.375D0, 741 * D3P875 = 3.875D0, DP25 = 0.25D0) 742 DATA C1,C3,C4,C5,C6 /1.D-15,1.D-20,1.D-14,1.D-9,1.D-5/ 743Cwas: DATA C1,C2,C3,C4,C5,C6/1.D-12,1.D-12,1.D-20,1.D-14,1.D-9,1.D-5/ 744 745 IF (NB.LE.1 .OR. NMAX.LE.0) RETURN 746 747 CALL QENTER('JACO_THR') 748c CALL GETTIM(TSTRT, WSTRT) 749 750 NNB = (NB*NB+NB)/2 751 752C 753C 921030-hjaaj: SD = C1 now 754C 900727-hjaaj: 755C SD calculation was done in every Jacobi iteration. 756C Now the largest absolute element in F is found once and 757C SD = 1.05D0 758C DO J = 1,NNB 759C SD = MAX(SD, ABS(F(J)) ) 760C END DO 761C SD=MAX(C1,C2*SD) 762 SD = max(C1, THR_conv) 763 764 JBIG(1:NB) = 0 765 BIG(1:NB) = D0 766 DO I = 1,NB 767 J = MIN(I-1,NMAX) 768 IF (J .GT. 0) THEN 769 II = (I*I-I)/2 770 JBIGI = 0 771 ABIGI = SD 772 DO K=1,J 773 IF (ABIGI .GE. ABS(F(II+K))) CYCLE 774 ABIGI = ABS(F(II+K)) 775 JBIGI = K 776 END DO 777 IF (JBIGI .GT. 0) THEN 778 JBIG(I) = JBIGI 779 BIG(I) = F(II+JBIGI) 780 END IF 781 END IF 782 END DO 783 784 MXITJA = 50*NNB 785 786 ITJACO = 0 787 200 ITJACO = ITJACO + 1 788 IF (ITJACO .GT. MXITJA) THEN 789 CALL QUIT('ERROR: JACO_THR did not converge ...') 790 END IF 791 792 T = D0 793 DO I = 2,NB 794 IF (T .GE. ABS(BIG(I))) CYCLE 795 T = ABS(BIG(I)) 796 IB = I 797 END DO 798 799 IF (T .LT. SD) GO TO 8000 ! finished 800 801 IA = JBIG(IB) 802 IAA = IA*(IA-1)/2 803 IBB = IB*(IB-1)/2 804 DIF = F(IAA+IA) - F(IBB+IB) 805 IF( ABS(DIF) .LE. C3) THEN 806 CX = ROOT2 807 SX = ROOT2 808 ELSE 809 T2X2 = BIG(IB)/DIF 810 T2X25 = T2X2*T2X2 811 IF(T2X25 .GT. C4) GO TO 240 812 CX = D1 813 SX = T2X2 814 GO TO 270 815 240 IF(T2X25 .GT. C5) GO TO 250 816 CX = D1 - DP5*T2X25 817 SX = T2X2*(D1 - D1P5*T2X25) 818 GO TO 270 819 250 IF(T2X25 .GT. C6) GO TO 260 820 CX = D1 + T2X25*(T2X25*D1P375 - DP5 ) 821 SX = T2X2*(D1 + T2X25*(T2X25*D3P875 - D1P5)) 822 GO TO 270 823 260 T = DP25 / SQRT(DP25 + T2X25) 824 CX = SQRT(DP5 + T) 825 SX = SIGN( SQRT(DP5 - T),T2X2 ) 826 270 CONTINUE 827 END IF 828 829 DO 275 IR=1,IA 830 T = F(IAA+IR)*SX 831 F(IAA+IR)= F(IAA+IR)*CX+F(IBB+IR)*SX 832 F(IBB+IR)=-T +F(IBB+IR)*CX 833 275 CONTINUE 834 IEAA=IAA+IA 835 IEAB=IBB+IA 836 TT =F(IEAB) 837 F(IEAB)=BIG(IB) 838 IF (JBIG(IA) .EQ. 0) THEN 839 IRST = IA + 1 840 IEAR = IEAA + IA 841 IEBR = IEAB + 1 842 ELSE 843 IRST = IA 844 IEAR = IEAA 845 IEBR = IEAB 846 END IF 847 DO 390 IR = IRST,NB 848 IF (IR .EQ. IA) GO TO 360 849C ... we have checked above that JBIG(IA) .ne. 0 850 T = F(IEAR)*SX 851 F(IEAR)= F(IEAR)*CX+F(IEBR)*SX 852 F(IEBR)=-T +F(IEBR)*CX 853 T =F(IEAR) 854 IT =IA 855 IF(IR-IB) 340,310,320 856 310 F(IEAA)=F(IEAA)*CX+F(IEAB)*SX 857C 921030+hjaaj: zero f(ab) to avoid round-off errors 858C F(IEAB)= TT*CX+F(IEBR)*SX 859 F(IEAB)= D0 860 F(IEBR)= -TT*SX+F(IEBR)*CX 861 GO TO 360 862 320 IF(ABS(T) .GE. ABS(F(IEBR))) GO TO 340 863 T =F(IEBR) 864 IT =IB 865 340 IF(ABS(T) .LT. ABS(BIG(IR))) GO TO 350 866 BIG(IR) = T 867 JBIG(IR) = IT 868 GO TO 380 869 350 IF(IA .NE. JBIG(IR) .AND. IB .NE. JBIG(IR)) GO TO 380 870 360 K= IEAR - IA 871 JBIGI = 0 872 IR1=MIN (IR-1,NMAX) 873 IF (IR1 .GT. 0) THEN 874 ABIGI = SD 875 DO I=1,IR1 876 IF(ABIGI .GE. ABS(F(K+I))) CYCLE 877 ABIGI = ABS(F(K+I)) 878 JBIGI =I 879 END DO 880 END IF 881 IF (JBIGI .GT. 0) THEN 882 JBIG(IR) = JBIGI 883 BIG(IR) = F(K+JBIGI) 884 ELSE 885 JBIG(IR) = 0 886 BIG(IR) = D0 887 END IF 888 380 CONTINUE 889 IEAR = IEAR + IR 890 IF (IR .GE. IB) THEN 891 IEBR = IEBR + IR 892 ELSE 893 IEBR = IEBR + 1 894 END IF 895 390 CONTINUE 896C 897 DO I=1,NROWV 898 VIIA = V(I,IA) 899 VIIB = V(I,IB) 900 V(I,IA) = VIIA*CX + VIIB*SX 901 V(I,IB) = -VIIA*SX + VIIB*CX 902 END DO 903 GO TO 200 904 905 8000 CONTINUE 906c CALL GETTIM(TEND, WEND) 907c WRITE(LUPRI,'(/A,4I10/A,2F20.2)') 908c & 'JACO_THR -- ITJACO, NB,NMAX,NROWV :',ITJACO, NB,NMAX,NROWV, 909c & 'JACO_THR -- CPUTIME, WALLTIME :',TEND-TSTRT, WEND-WSTRT 910 CALL QEXIT('JACO_THR') 911 RETURN 912 END 913C /* Deck norm */ 914 SUBROUTINE NORM(S,VC,N,M,W,THNORM,IRETUR) 915C 916C revised 14-May-1985 hjaaj (call MPAPV instead of CNTRC) 917C revised May 2000 hjaaj (two rounds if small norm) 918C 919C COMPUTES SCHMIDT-ORTHONORMALIZED SET OF VECTORS 920C CALLING SEQUENCE PARAMETERS ARE AS FOLLOWS 921C S METRIC MATRIX STORED AS LOWER TRIANGLE (R*8) 922C VC LOCATION OF ORIGINAL NON-ORTHONORMAL VECTORS (R*8) 923C FINAL ORTHONORMALIZED VECTORS REPLACE ORIGINAL SET 924C N DIMENSION OF BASIS SET (I*4) 925C M NUMBER OF VECTORS TO BE ORTHONORMALIZED (I*4) 926C W TEMPORARY WORKING AREA OF 2*N WORDS (R*8) 927C RETURNS 928C NORMAL RETURN ORTHONORMALIZED SET OBTAINED 929C RETURN 1 INITIAL VECTORS AT VC LINEARLY 930C DEPENDENT WITHIN THRES HOLD (THNORM) 931C AUXILIARY ENTRY 932C 933#include "implicit.h" 934#include "priunit.h" 935 DIMENSION S(*), VC(N,M), W(*) 936 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, THRRND = 0.9D0 ) 937 IRETUR=0 938C 939C N = 1 special case 940C 941 IF (N .EQ. 1) THEN 942 IF (M .EQ. 1) THEN 943 IF (VC(1,1)*VC(1,1) .LT. THNORM) THEN 944 IRETUR=-1 945 ELSE 946 VC(1,1) = D1/SQRT(S(1)) 947 END IF 948 END IF 949 RETURN 950 END IF 951C 952C BEGIN OUTERMOST LOOP OVER TRIAL VECTOR SET 953C 954 DO 20 I=1,M 955 ITURN = 0 956 1 ITURN = ITURN + 1 957 IROUND = 0 958C 959 CALL MPAPV(N,S,VC(1,I),W) 960 TNORM = DDOT(N,VC(1,I),1,W(1),1) 961 IF (TNORM .LT. THNORM) THEN 962Chj ... zero vector on input 963 IRETUR = -I 964 RETURN 965 END IF 966Chj may2000: normalize input vector 967C (we ignore round-off errors as it is renormalized later) 968 TNORM = D1 / SQRT(TNORM) 969 CALL DSCAL(N,TNORM,VC(1,I),1) 970 CALL MPAPV(N,S,VC(1,I),W) 971 TNORM = DDOT(N,VC(1,I),1,W(1),1) 972C 973C BEGIN COEFFICIENTS AND NORMALIZATION LOOP 974C 975 DO 5 J=1,I-1 976 T = DDOT(N,VC(1,J),1,W(1),1) 977 TNORM = TNORM - T*T 978 5 W(N+J) = -T 979 IF (TNORM .LT. THNORM) THEN 980Chj ... zero vector after orthogonalization 981 IRETUR = I 982 RETURN 983 END IF 984 IF (TNORM .LT. THRRND) IROUND = IROUND + 1 985Chj ... experiments have shown that TNORM as big as 986C 0.25 can give a normalization error of 1.0D-7 ! 987 TNORM = D1/SQRT(TNORM) 988 DO J=1,I-1 989 W(N+J) = W(N+J)*TNORM 990 END DO 991 W(N+I) = TNORM 992C 993C REPLACE VC(*,I) 994C 995 CALL DGEMM('N','N',N,1,I,1.D0, 996 & VC,N, 997 & W(N+1),I,0.D0, 998 & W,N) 999 CALL DCOPY(N,W,1,VC(1,I),1) 1000C 1001 IF (ITURN .EQ. 1 .AND. IROUND .GT. 0) THEN 1002 IF (IPRSTAT .GT. 0) THEN 1003 WRITE (LUSTAT,*) 'Info: second round in NORM, I=',I 1004 CALL QDUMP(LUSTAT) 1005 END IF 1006 GO TO 1 1007 END IF 1008 IF (IROUND.GT.0) CALL QUIT('NORM: round-off errors, see code') 1009Chj ... this ought never to happen ... 1010 20 CONTINUE 1011 RETURN 1012 END 1013C /* Deck readi */ 1014 SUBROUTINE READI (IT,N,INTX) 1015C 1016C (30-Jan-1984) hjaaj 1017C 1018 INTEGER IT, N, INTX(N) 1019 IF (N .GT. 0) THEN 1020 READ (IT,END = 10) INTX 1021 ELSE 1022 READ (IT,END = 10) 1023 END IF 1024 RETURN 1025 10 CONTINUE 1026 INTX(N)=-1 1027 RETURN 1028 END 1029C /* Deck readi4 */ 1030 SUBROUTINE READI4(IT,N,INTX) 1031C 1032C (Jul-2014) hjaaj, based on READI 1033C Purpose: when INTX is always INTEGER*4 1034C 1035 INTEGER IT, N 1036 INTEGER*4 INTX(N) 1037 IF (N .GT. 0) THEN 1038 READ (IT,END = 10) INTX 1039 ELSE 1040 READ (IT,END = 10) 1041 END IF 1042 RETURN 1043 10 CONTINUE 1044 INTX(N)=-1 1045 RETURN 1046 END 1047C /* Deck readi8 */ 1048 SUBROUTINE READI8(IT,N,INTX) 1049C 1050C (Jul-2014) hjaaj, based on READI 1051C Purpose: when INTX is always INTEGER*8 1052C 1053 INTEGER IT, N 1054 INTEGER*8 INTX(N) 1055 IF (N .GT. 0) THEN 1056 READ (IT,END = 10) INTX 1057 ELSE 1058 READ (IT,END = 10) 1059 END IF 1060 RETURN 1061 10 CONTINUE 1062 INTX(N)=-1 1063 RETURN 1064 END 1065C /* Deck readdi */ 1066 SUBROUTINE READDI(IT,IU,N,IX) 1067 DIMENSION IX(N) 1068 READ(IT,REC=IU) IX 1069 RETURN 1070 END 1071C /* Deck readt */ 1072 SUBROUTINE READT (IT,N,X) 1073#include "implicit.h" 1074#include "priunit.h" 1075 CHARACTER*120 FNAME 1076 DIMENSION X(N) 1077 IF (IT .LE. 0) GOTO 30 1078 READ (IT,IOSTAT=IERR,END=10,ERR=20) X 1079 RETURN 1080 10 CONTINUE 1081 INQUIRE(UNIT=IT,NAME=FNAME) 1082 WRITE (LUPRI,*) 'READT: END reading file ',FNAME 1083 WRITE (LUPRI,*) 'UNIT ',IT, ' record length ',N 1084 WRITE ( 0 ,*) 'READT: END reading file ',FNAME 1085 WRITE ( 0 ,*) 'UNIT ',IT, ' record length ',N 1086 CALL QUIT('READT: END reading file') 1087 20 CONTINUE 1088 INQUIRE(UNIT=IT,NAME=FNAME) 1089 WRITE (LUPRI,*) 'READT: ERROR reading file ',FNAME 1090 WRITE (LUPRI,*) 'UNIT ',IT, ' record length ',N,' error code ',IERR 1091 WRITE ( 0 ,*) 'READT: ERROR reading file ',FNAME 1092 WRITE ( 0 ,*) 'UNIT ',IT, ' record length ',N,' error code ',IERR 1093 CALL QUIT('READT: Error reading file') 1094 30 CONTINUE 1095 WRITE (LUPRI,*) 'READT ERRROR: non-positive file unit number ',IT 1096 CALL QUIT('READT ERRROR: non-positive file unit number') 1097 END 1098C /* Deck writi */ 1099 SUBROUTINE WRITI (IT,N,INTX) 1100C 1101C (30-Jan-1984) hjaaj 1102C 1103 INTEGER IT, N, INTX(N) 1104 WRITE (IT) INTX 1105 RETURN 1106 END 1107C /* Deck writi4 */ 1108 SUBROUTINE WRITI4(IT,N,INTX) 1109C 1110C (May-2016) hjaaj, based on WRITI 1111C Purpose: when INTX is always INTEGER*4 1112C 1113 INTEGER IT, N 1114 INTEGER*4 INTX(N) 1115 WRITE (IT) INTX 1116 RETURN 1117 END 1118C /* Deck writi8 */ 1119 SUBROUTINE WRITI8(IT,N,INTX) 1120C 1121C (May-2016) hjaaj, based on WRITI 1122C Purpose: when INTX is always INTEGER*8 1123C 1124 INTEGER IT, N 1125 INTEGER*8 INTX(N) 1126 WRITE (IT) INTX 1127 RETURN 1128 END 1129#if defined (VAR_SPLITFILES) 1130C /* Deck writsi */ 1131 SUBROUTINE WRITSI(IT,N,INTX,JBUF) 1132C 1133C K.Ruud, April 13 2000 1134C 1135#include "implicit.h" 1136#include "2gbdef.h" 1137#include "priunit.h" 1138 DIMENSION INTX(N) 1139 CHARACTER*80 FNNAME, FNNM2 1140#include "inftap.h" 1141#include "chrnos.h" 1142C 1143 IF ((JBUF + N) .GT. I2GB) THEN 1144C 1145C Ooops, file will be overfull, need to open a new one...... 1146C 1147 INQUIRE(UNIT=IT,NAME=FNNAME) 1148 LN = 1 1149 10 CONTINUE 1150 IF (FNNAME(LN:LN) .NE. ' ') THEN 1151 LN = LN + 1 1152 GOTO 10 1153 END IF 1154 LN = LN - 1 1155 CALL GPCLOSE(IT,'KEEP') 1156 I = LN - 1 1157 IF (FNNAME(I:I) .NE. '-') THEN 1158 FNNM2 = FNNAME(1:LN)//'-0' 1159 LN = LN + 2 1160 ELSE 1161 READ(FNNAME(LN:),'(I1)') INUM 1162 INUM = INUM + 1 1163 IF (INUM .GT. 9) THEN 1164 WRITE (LUPRI,'(/A)') ' DALTON needs to split a file '// 1165 & ' more than 11 times.', 1166 & ' This is currently not supported' 1167 CALL QUIT('Too many splittings of a file') 1168 END IF 1169 FNNM2 = FNNAME(1:I)//CHRNOS(INUM) 1170 END IF 1171 CALL GPOPEN(IT,FNNM2(1:LN),'UNKNOWN',' ',' ',IDUMMY, 1172 & .FALSE.) 1173 JBUF = 0 1174 END IF 1175 WRITE (IT) INTX 1176 JBUF = JBUF + N 1177 RETURN 1178 END 1179 SUBROUTINE WRITST(IT,N,X,JBUF) 1180C 1181C K.Ruud, April 13 2000 1182C 1183#include "implicit.h" 1184#include "2gbdef.h" 1185#include "priunit.h" 1186 DIMENSION X(N) 1187 CHARACTER*80 FNNAME, FNNM2 1188#include "inftap.h" 1189#include "chrnos.h" 1190C 1191 WRITE (IT) N 1192 IF ((JBUF + N + 1) .GT. I2GB) THEN 1193C 1194C Ooops, file will be overfull, need to open a new one...... 1195C 1196 INQUIRE(UNIT=IT,NAME=FNNAME) 1197 LN = 1 1198 10 CONTINUE 1199 IF (FNNAME(LN:LN) .NE. ' ') THEN 1200 LN = LN + 1 1201 GOTO 10 1202 END IF 1203 LN = LN - 1 1204 CALL GPCLOSE(IT,'KEEP') 1205 I = LN - 1 1206 IF (FNNAME(I:I) .NE. '-') THEN 1207 FNNM2 = FNNAME(1:LN)//'-0' 1208 LN = LN + 2 1209 ELSE 1210 READ(FNNAME(LN:),'(I1)') INUM 1211 INUM = INUM + 1 1212 IF (INUM .GT. 9) THEN 1213 WRITE (LUPRI,'(/A)') ' DALTON needs to split a file '// 1214 & ' more than 11 times.', 1215 & ' This is currently not supported' 1216 CALL QUIT('Too many splittings of a file') 1217 END IF 1218 FNNM2 = FNNAME(1:I)//CHRNOS(INUM) 1219 END IF 1220 CALL GPOPEN(IT,FNNM2(1:LN),'UNKNOWN',' ',' ',IDUMMY, 1221 & .FALSE.) 1222 JBUF = 0 1223 END IF 1224 WRITE (IT) X 1225 JBUF = JBUF + N + 1 1226 RETURN 1227 END 1228C /* Deck readsi */ 1229 SUBROUTINE READSI(IT,N,INTX,JBUF) 1230C 1231C K.Ruud, April 13 2000 1232C 1233#include "implicit.h" 1234#include "2gbdef.h" 1235#include "priunit.h" 1236 DIMENSION INTX(N) 1237 CHARACTER*80 FNNAME, FNNM2 1238#include "inftap.h" 1239#include "chrnos.h" 1240C 1241 IF ((JBUF + N) .GT. I2GB) THEN 1242C 1243C Ooops, file will be overfull, need to open a new one...... 1244C 1245 INQUIRE(UNIT=IT,NAME=FNNAME) 1246 LN = 1 1247 10 CONTINUE 1248 IF (FNNAME(LN:LN) .NE. ' ') THEN 1249 LN = LN + 1 1250 GOTO 10 1251 END IF 1252 LN = LN - 1 1253 CALL GPCLOSE(IT,'KEEP') 1254 I = LN - 1 1255 IF (FNNAME(I:I) .NE. '-') THEN 1256 FNNM2 = FNNAME(1:LN)//'-0' 1257 LN = LN + 2 1258 ELSE 1259 READ(FNNAME(LN:),'(I1)') INUM 1260 INUM = INUM + 1 1261 IF (INUM .GT. 9) THEN 1262 WRITE (LUPRI,'(/A)') ' DALTON wants to read from a '// 1263 & ' file split more than 11 times.', 1264 & ' This is currently not supported' 1265 CALL QUIT('Too many splittings of a file') 1266 END IF 1267 FNNM2 = FNNAME(1:I)//CHRNOS(INUM) 1268 END IF 1269 CALL GPOPEN(IT,FNNM2(1:LN),'UNKNOWN',' ',' ',IDUMMY, 1270 & .FALSE.) 1271 REWIND (IT) 1272 JBUF = 0 1273 END IF 1274 READ (IT, END = 30) INTX 1275 JBUF = JBUF + N 1276 RETURN 1277 30 INTX(N) = -1 1278 JBUF = JBUF + 1 1279 RETURN 1280 END 1281 SUBROUTINE READST(IT,N,X,JBUF,READ) 1282C 1283C K.Ruud, Dec 08 2000 1284C 1285#include "implicit.h" 1286#include "2gbdef.h" 1287#include "priunit.h" 1288 DIMENSION X(*) 1289 LOGICAL READ 1290 CHARACTER*80 FNNAME, FNNM2 1291#include "inftap.h" 1292#include "chrnos.h" 1293C 1294 READ (IT, END=30) N 1295 IF ((JBUF + N + 1) .GT. I2GB) THEN 1296C 1297C Ooops, file will be overfull, need to open a new one...... 1298C 1299 INQUIRE(UNIT=IT,NAME=FNNAME) 1300 LN = 1 1301 10 CONTINUE 1302 IF (FNNAME(LN:LN) .NE. ' ') THEN 1303 LN = LN + 1 1304 GOTO 10 1305 END IF 1306 LN = LN - 1 1307 CALL GPCLOSE(IT,'KEEP') 1308 I = LN - 1 1309 IF (FNNAME(I:I) .NE. '-') THEN 1310 FNNM2 = FNNAME(1:LN)//'-0' 1311 LN = LN + 2 1312 ELSE 1313 READ(FNNAME(LN:),'(I1)') INUM 1314 INUM = INUM + 1 1315 IF (INUM .GT. 9) THEN 1316 WRITE (LUPRI,'(/A)') ' DALTON wants to read from a '// 1317 & ' file split more than 11 times.', 1318 & ' This is currently not supported' 1319 CALL QUIT('Too many splittings of a file') 1320 END IF 1321 FNNM2 = FNNAME(1:I)//CHRNOS(INUM) 1322 END IF 1323 CALL GPOPEN(IT,FNNM2(1:LN),'UNKNOWN',' ',' ',IDUMMY, 1324 & .FALSE.) 1325 REWIND (IT) 1326 JBUF = 0 1327 END IF 1328 IF (READ) THEN 1329 READ (IT, END = 30) (X(I),I = 1, N) 1330 ELSE 1331 READ (IT, END = 30) 1332 END IF 1333 JBUF = JBUF + N + 1 1334 RETURN 1335 30 N = -1 1336 JBUF = JBUF + 1 1337 RETURN 1338 END 1339#endif 1340C /* Deck writdi */ 1341 SUBROUTINE WRITDI(IT,IU,N,IX) 1342 DIMENSION IX(N) 1343 WRITE(IT,REC=IU) IX 1344 RETURN 1345 END 1346C /* Deck writt */ 1347 SUBROUTINE WRITT (IT,N,X) 1348#include "implicit.h" 1349 DIMENSION X(N) 1350 WRITE (IT) X 1351 RETURN 1352C 1353 END 1354C /* Deck mollab */ 1355 SUBROUTINE MOLLAB(A,LU,LUERR) 1356C 1357C 16-Jun-1986 hjaaj 1358C (as SEARCH but CHARACTER*8 instead of REAL*8 labels) 1359C 1360C Purpose: 1361C Search for MOLECULE labels on file LU 1362C 1363 CHARACTER*8 A, B(4), C 1364#if defined (VAR_SPLITFILES) 1365#include "dummy.h" 1366#endif 1367 CHARACTER FNNAME*80 1368 DATA C/'********'/ 1369 IRDERR = 0 1370#if defined (VAR_SPLITFILES) 1371 INQUIRE (UNIT=LU,NAME=FNNAME) 1372 LN = 1 1373 10 CONTINUE 1374 IF (FNNAME(LN:LN) .EQ. '-') THEN 1375 LN = LN - 1 1376 LUBKP = LU 1377 CALL GPCLOSE(LU,'KEEP') 1378 LU = LUBKP 1379 CALL GPOPEN(LU,FNNAME(1:LN),'OLD','SEQUENTIAL','UNFORMATTED', 1380 & IDUMMY,.FALSE.) 1381 INQUIRE (UNIT=LU,NAME=FNNAME) 1382 REWIND (LU) 1383 GOTO 1 1384 ELSE IF (FNNAME(LN:LN) .EQ. ' ') THEN 1385 GOTO 1 1386 END IF 1387 LN = LN + 1 1388 GOTO 10 1389#endif 1390 1 READ (LU,END=3,ERR=6,IOSTAT=IOSVAL) B 1391 IRDERR = 0 1392 IF (B(1).NE.C) GO TO 1 1393 IF (B(4).NE.A) GO TO 1 1394 IF (LUERR.LT.0) LUERR = 0 1395 RETURN 1396C 1397 3 IF (LUERR.LT.0) THEN 1398#if defined (VAR_MFDS) 1399C 880916-hjaaj -- attempt to fix an IBM problem 1400C IBM shifts to new file after END= branch (e.g. FTxxF002 instead 1401C of FTxxF001), backspace makes LU ready for append. 1402C 940510-hjaaj: same change for Cray multifile datasets 1403 BACKSPACE LU 1404#endif 1405 LUERR = -1 1406 RETURN 1407 ELSE 1408 INQUIRE (UNIT=LU,NAME=FNNAME) 1409 WRITE(LUERR,4)A,LU,FNNAME 1410 CALL QTRACE(LUERR) 1411 REWIND (LU) 1412 CALL DMPLAB(LU,LUERR) 1413 CALL QUIT('ERROR (MOLLAB) MOLECULE label not found on file') 1414 END IF 1415 4 FORMAT(/' *** ERROR (MOLLAB), MOLECULE label ',A8, 1416 * ' not found on unit',I4/' File name: ',A) 1417C 1418 6 IRDERR = IRDERR + 1 1419 IF (IRDERR .LT. 10) GO TO 1 1420 IF (LUERR.LT.0) THEN 1421 LUERR = -2 1422 RETURN 1423 ELSE 1424 WRITE (LUERR,7) LU,A,IOSVAL 1425 CALL QTRACE(LUERR) 1426 CALL QUIT('ERROR (MOLLAB) error reading file') 1427 END IF 1428 7 FORMAT(/' *** ERROR (MOLLAB), error reading unit',I4, 1429 * /T22,'when searching for label ',A8, 1430 * /T22,'IOSTAT value :',I7) 1431 END 1432C /* Deck fndlab */ 1433 LOGICAL FUNCTION FNDLAB(A,LU) 1434C 1435C 26-May-1985 hjaaj -- logical function version of SEARCH 1436C 16-Jun-1986 hjaaj -- changed to CHARACTER*8 from REAL*8 1437C 1438#include "priunit.h" 1439 CHARACTER*8 A, B(4), C 1440 DATA C/'********'/ 1441 IRDERR = 0 1442 1 READ(LU,END=3,ERR=6)B 1443 IRDERR = 0 1444 IF (B(1).NE.C) GO TO 1 1445 IF (B(4).NE.A) GO TO 1 1446 FNDLAB = .TRUE. 1447 GO TO 10 1448C 1449 6 IRDERR = IRDERR + 1 1450 IF (IRDERR .LT. 5) GO TO 1 1451 GO TO 8 1452 3 CONTINUE 1453#if defined (VAR_MFDS) 1454C 880916-hjaaj -- attempt to fix an IBM problem 1455C IBM shifts to new file after END= branch (e.g. FTxxF002 instead 1456C of FTxxF001), backspace makes LU ready for append. 1457C 940510-hjaaj: same change for Cray multifile datasets 1458 BACKSPACE (LU,ERR=7) 1459C ... aug07-hjaaj: new ERR branch to avoid problems with backspace on empty files. 1460 GO TO 8 1461 7 write (LUPRI,*) 1462 & 'FNDLAB WARNING: ERR branch for backspace on LU ',LU 1463 call qdump(LUPRI) 1464#endif 1465C 1466 8 FNDLAB = .FALSE. 1467C 1468 10 RETURN 1469 END 1470C /* Deck mollb2 */ 1471 SUBROUTINE MOLLB2(SRCLBL,RTNLBL,LU,LUERR) 1472C 1473C 28-Jun-1986 hjaaj 1474C (as MOLLAB, but returns two middle labels in RTNLBL(2)) 1475C 1476C Purpose: 1477C Search for MOLECULE labels on file LU 1478C 1479 CHARACTER*8 SRCLBL, RTNLBL(2), B(4), STAR8 1480 PARAMETER (STAR8 = '********') 1481C 1482 IRDERR = 0 1483 1 READ (LU,END=3,ERR=6,IOSTAT=IOSVAL) B 1484 IRDERR = 0 1485 IF (B(1).NE.STAR8) GO TO 1 1486 IF (B(4).NE.SRCLBL) GO TO 1 1487C 1488 RTNLBL(1) = B(2) 1489 RTNLBL(2) = B(3) 1490 IF (LUERR.LT.0) LUERR = 0 1491 RETURN 1492C 1493 3 IF (LUERR.LT.0) THEN 1494#if defined (VAR_MFDS) 1495C 880916-hjaaj -- attempt to fix an IBM problem 1496C IBM shifts to new file after END= branch (e.g. FTxxF002 instead 1497C of FTxxF001), backspace makes LU ready for append. 1498C 940510-hjaaj: same change for Cray multifile datasets 1499 BACKSPACE LU 1500#endif 1501 LUERR = -1 1502 RETURN 1503 ELSE 1504 WRITE (LUERR,4) SRCLBL,LU 1505 CALL QTRACE(LUERR) 1506 CALL QUIT('ERROR (MOLLB2) MOLECULE label not found on file') 1507 END IF 1508 4 FORMAT(/' *** ERROR (MOLLB2), MOLECULE label ',A8, 1509 * ' not found on unit',I4) 1510C 1511 6 IRDERR = IRDERR + 1 1512 IF (IRDERR .LT. 5) GO TO 1 1513 IF (LUERR.LT.0) THEN 1514 LUERR = -2 1515 RETURN 1516 ELSE 1517 WRITE (LUERR,7) LU,SRCLBL,IOSVAL 1518 CALL QTRACE(LUERR) 1519 CALL QUIT('ERROR (MOLLB2) error reading file') 1520 END IF 1521 7 FORMAT(/' *** ERROR (MOLLB2), error reading unit',I4, 1522 * /T22,'when searching for label ',A8, 1523 * /T22,'IOSTAT value :',I7) 1524 END 1525C /* Deck fndlb2 */ 1526 LOGICAL FUNCTION FNDLB2(SRCLBL,RTNLBL,LU) 1527C 1528C 5-Aug-1986 hjaaj 1529C (as FNDLAB, but returns two middle labels in RTNLBL(2)) 1530C 1531 CHARACTER*8 SRCLBL, RTNLBL(2), B(4), STAR8 1532 PARAMETER (STAR8 = '********') 1533 IRDERR = 0 1534 1 READ (LU,END=3,ERR=6) B 1535 IRDERR = 0 1536 IF (B(1).NE.STAR8) GO TO 1 1537 IF (B(4).NE.SRCLBL) GO TO 1 1538 FNDLB2 = .TRUE. 1539 RTNLBL(1) = B(2) 1540 RTNLBL(2) = B(3) 1541 GO TO 10 1542C 1543 6 IRDERR = IRDERR + 1 1544 IF (IRDERR .LT. 5) GO TO 1 1545 GO TO 8 1546 3 CONTINUE 1547#if defined (VAR_MFDS) 1548C 880916-hjaaj -- attempt to fix an IBM problem 1549C IBM shifts to new file after END= branch (e.g. FTxxF002 instead 1550C of FTxxF001), backspace makes LU ready for append. 1551C 940510-hjaaj: same change for Cray multifile datasets 1552 BACKSPACE LU 1553#endif 1554 8 FNDLB2 = .FALSE. 1555C 1556 10 RETURN 1557 END 1558C /* Deck fndlb3 */ 1559 SUBROUTINE FNDLB3(SRCLBL,IVALUE,LU) 1560C 1561C 06-Jun-2000 ALig 1562C (as FNDLB2, but the IVALUE returns the value of the 1563C symmetry representation of the operator SRCLBL and 0 if the 1564C label does not exist) 1565C 1566#include "implicit.h" 1567 CHARACTER*8 SRCLBL, B(4), STAR8 1568 PARAMETER (STAR8 = '********') 1569 IRDERR = 0 1570 REWIND(LU) 1571 1 READ (LU,END=3,ERR=6) B 1572 IRDERR = 0 1573 IF (B(1).NE.STAR8) GO TO 1 1574 IF (B(4).NE.SRCLBL) GO TO 1 1575 IVALUE = (ICHAR(B(2)(2:2))-ICHAR('0')) 1576 GO TO 10 1577C 1578 6 IRDERR = IRDERR + 1 1579 IF (IRDERR .LT. 5) GO TO 1 1580 GO TO 8 1581 3 CONTINUE 1582#if defined (VAR_MFDS) 1583C 880916-hjaaj -- attempt to fix an IBM problem 1584C IBM shifts to new file after END= branch (e.g. FTxxF002 instead 1585C of FTxxF001), backspace makes LU ready for append. 1586C 940510-hjaaj: same change for Cray multifile datasets 1587 BACKSPACE LU 1588#endif 1589 8 IVALUE = 0 1590C 1591 10 CONTINUE 1592 END 1593C 1594C /* Deck nxtlab */ 1595 LOGICAL FUNCTION NXTLAB(SRCLBL, RTNLBL, LU) 1596C 1597C 3-Nov-1986 hjaaj 1598C (find and return next MOLECULE label on LU, 1599C NXTLAB false if no label found) 1600C 1601 CHARACTER*8 SRCLBL, RTNLBL(2), B(4), STAR8 1602 PARAMETER ( STAR8 = '********' ) 1603 IRDERR = 0 1604 1 READ(LU,END=3,ERR=6) B 1605 IRDERR = 0 1606 IF (B(1) .NE. STAR8) GO TO 1 1607 NXTLAB = .TRUE. 1608 SRCLBL = B(4) 1609 RTNLBL(1) = B(2) 1610 RTNLBL(2) = B(3) 1611 GO TO 10 1612C 1613 6 IRDERR = IRDERR + 1 1614 IF (IRDERR .LT. 5) GO TO 1 1615 GO TO 8 1616 3 CONTINUE 1617#if defined (VAR_MFDS) 1618C 880916-hjaaj -- attempt to fix an IBM problem 1619C IBM shifts to new file after END= branch (e.g. FTxxF002 instead 1620C of FTxxF001), backspace makes LU ready for append. 1621C 940510-hjaaj: same change for Cray multifile datasets 1622 BACKSPACE LU 1623#endif 1624 8 NXTLAB = .FALSE. 1625C 1626 10 RETURN 1627 END 1628C /* Deck dmplab */ 1629 SUBROUTINE DMPLAB(LU,LUPRI) 1630C 1631C 27-Mar-1987 hjaaj -- dump remaining labels on file LU 1632C 1633 CHARACTER*8 B(4), C 1634 PARAMETER ( C = '********' ) 1635C 1636 WRITE (LUPRI, '(//A,I5)') ' --- DUMP OF LABELS ON UNIT',LU 1637 IRDERR = 0 1638 IREC = 0 1639 1 READ (LU,END=3,ERR=6,IOSTAT=IOSVAL) B 1640 IRDERR = 0 1641 IREC = IREC + 1 1642 IF (B(1).EQ.C) THEN 1643 WRITE (LUPRI, '(A,I6,4(2X,A8))') ' Rec. no.',IREC,B 1644 END IF 1645 GO TO 1 1646C 1647 6 CONTINUE 1648 IRDERR = IRDERR + 1 1649 IREC = IREC + 1 1650 WRITE (LUPRI, '(/A,I6,A,I7)') 1651 & ' ERROR reading rec. no.',IREC,'; IOSTAT value',IOSVAL 1652 IF (IRDERR .LT. 5) GO TO 1 1653 WRITE (LUPRI, '(/A)') 1654 & ' ERROR exit from DMPLAB: 5 consecutive read errors ---' 1655 3 CONTINUE 1656 REWIND (LU) 1657C 1658 WRITE (LUPRI, '(/I10,A)') IREC,' records read from file.' 1659 RETURN 1660 END 1661C /* Deck newlab */ 1662 SUBROUTINE NEWLAB(LABEL,LU,LUERR) 1663C 1664C 29-Sep-1988 Hans Joergen Aa. Jensen 1665C 1666C Write new MOLECULE-type label to LU 1667C 1668 CHARACTER*8 LABEL, LTIME, LDATE, LSTARS 1669 DATA LSTARS /'********'/ 1670 CALL GETDAT(LDATE,LTIME) 1671 WRITE (LU,ERR=1000,IOSTAT=IOSVAL) LSTARS,LDATE,LTIME,LABEL 1672 RETURN 1673C 1674 1000 IF (LUERR .LT. 0) THEN 1675 LUERR = -2 1676 RETURN 1677 ELSE 1678 WRITE (LUERR,'(/3A,I5/A,I7)') 1679 & ' NEWLAB: error writing label "',LABEL,'" to unit',LU, 1680 & ' IOSTAT value',IOSVAL 1681 CALL QTRACE(LUERR) 1682 CALL QUIT('NEWLAB: output error writing label') 1683 END IF 1684 END 1685C /* Deck newlb2 */ 1686 SUBROUTINE NEWLB2(LABEL,RTNLBL,LU,LUERR) 1687C 1688C 29-Sep-1988 Hans Joergen Aa. Jensen 1689C 1690C Write new MOLECULE-type label to LU 1691C 1692 CHARACTER*8 LABEL, RTNLBL(2), LSTARS 1693 DATA LSTARS /'********'/ 1694 WRITE (LU,ERR=1000,IOSTAT=IOSVAL) LSTARS,RTNLBL,LABEL 1695 RETURN 1696C 1697 1000 IF (LUERR .LT. 0) THEN 1698 LUERR = -2 1699 RETURN 1700 ELSE 1701 WRITE (LUERR,'(/3A,I5/A,I7)') 1702 & ' NEWLB2: error writing label "',LABEL,'" to unit',LU, 1703 & ' IOSTAT value',IOSVAL 1704 CALL QTRACE(LUERR) 1705 CALL QUIT('NEWLB2: output error writing label') 1706 END IF 1707 END 1708C /* Deck second */ 1709#if !defined (VAR_GFORTRAN) 1710 REAL*8 FUNCTION SECOND () 1711 REAL*4 TIME_CPU 1712 CALL CPU_TIME(TIME_CPU) 1713 SECOND = TIME_CPU 1714 RETURN 1715 END 1716#endif 1717C /* Deck sotmat */ 1718 SUBROUTINE SOTMAT(NMO,UMO,IFAIL) 1719C 1720C 16-Feb-1986 Hans Jorgen Aa. Jensen 1721C 1722C Purpose: 1723C 1724C Construct the orbital transformation matrix 1725C for transforming a CAS CI vector using a sequence 1726C of single orbital transformations as described by 1727C Per-Ake Malmquist. 1728C 1729C The matrix is 1730C 1731C C + C = (1 - L) + U(inv); 1732C L U 1733C 1734C where L and U constitute the LU decomposition of the 1735C orthogonal orbital transformation matrix UMO. 1736C 1737#include "implicit.h" 1738 DIMENSION UMO(NMO,NMO) 1739#include "priunit.h" 1740 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 ) 1741 PARAMETER ( THRESH = 1.D-4 ) 1742C 1743 CALL QENTER('SOTMAT') 1744C 1745C STEP 1: The LU decomposition 1746C ============================ 1747C 1748 DO 220 K = 1,NMO 1749 X = UMO(K,K) 1750 IF (ABS(X) .LE. THRESH) GO TO 960 1751 X = D1 / X 1752 DO 120 I = K+1,NMO 1753 UMO(I,K) = UMO(I,K) * X 1754 120 CONTINUE 1755 DO 200 J = (K+1),NMO 1756 Y = UMO(K,J) 1757 DO 180 I = K+1,NMO 1758 UMO(I,J) = UMO(I,J) - UMO(I,K)*Y 1759 180 CONTINUE 1760 200 CONTINUE 1761 220 CONTINUE 1762C 1763C STEP 2: U inverse 1764C ================= 1765C 1766 DET = D1 1767 DO 300 K = 1,NMO 1768 DET = DET * UMO(K,K) 1769 300 CONTINUE 1770#if defined (VAR_SOTMATTEST) 1771C 1772C 1773 WRITE (LUPRI,'(//A,1P,D10.2)') ' SOTMAT: UMO determinant =',DET 1774 WRITE (LUPRI,'(//A)') ' SOTMAT: LU matrix' 1775 CALL OUTPUT (UMO,1,NMO,1,NMO,NMO,NMO,1,6) 1776#endif 1777 IF (ABS(DET) .LE. THRESH) GO TO 980 1778C 1779 DO 400 K = 1,NMO 1780 UIDIAG = D1 / UMO(K,K) 1781 UMO(K,K) = UIDIAG 1782 DO 380 J = 1,(K-1) 1783 SUM = D0 1784 DO 360 I = J,(K-1) 1785 SUM = SUM - UMO(J,I)*UMO(I,K) 1786 360 CONTINUE 1787 UMO(J,K) = UIDIAG * SUM 1788 380 CONTINUE 1789 400 CONTINUE 1790C 1791C STEP 3: construct 1 - L 1792C ======================= 1793C 1794 DO 600 K = 1,NMO 1795 DO 580 J = (K+1),NMO 1796 UMO(J,K) = -UMO(J,K) 1797 580 CONTINUE 1798 600 CONTINUE 1799C 1800C NORMAL AND ERROR RETURNS 1801C ======================== 1802C 1803 IFAIL = 0 1804 GO TO 9999 1805C 1806 960 CONTINUE 1807 WRITE (LUPRI,'(///A,1P,D10.2/)') 1808 1 ' --- SOTMAT: DIAGONAL ELEMENT TOO SMALL:',X 1809 IFAIL = 1 1810 GO TO 9999 1811C 1812 980 CONTINUE 1813 WRITE (LUPRI,'(///2A,1P,D10.2/)') 1814 1 ' --- SOTMAT: U MATRIX TOO CLOSE TO SINGULARITY,', 1815 2 ' DETERMINANT =',DET 1816 IFAIL = 2 1817C 1818 9999 CALL QEXIT('SOTMAT') 1819 RETURN 1820C 1821C END OF SOTMAT. 1822C 1823 END 1824C /* Deck nofdia */ 1825 INTEGER FUNCTION NOFDIA(N,NDIM,AMAT,THREQL) 1826C 1827C 2-Jul-1992 hjaaj 1828C This function returns the number of off-diagonal elements 1829C with absolute value greater than THREQL. 1830C 1831#include "implicit.h" 1832 DIMENSION AMAT(NDIM,NDIM) 1833C 1834 CALL QENTER('NOFDIA') 1835 NELEM = 0 1836 DO 220 K = 1,N 1837 DO 210 I = 1,N 1838 IF (I .NE. K .AND. ABS(AMAT(I,K)) .GT. THREQL) THEN 1839 NELEM = NELEM + 1 1840 END IF 1841 210 CONTINUE 1842 220 CONTINUE 1843 NOFDIA = NELEM 1844 CALL QEXIT('NOFDIA') 1845 RETURN 1846 END 1847C /* Deck matsym */ 1848 INTEGER FUNCTION MATSYM(N,NDIM,AMAT,THRZER) 1849C 1850C 27-Apr-1999 hjaaj 1851C This function returns 1852C 1 if the matrix is symmetric to threshold THRZER 1853C 2 if the matrix is antisymmetric to threshold THRZER 1854C 3 if all elements are below THRZER 1855C 0 otherwise (the matrix is unsymmetric about the diagonal) 1856C 1857#include "implicit.h" 1858 DIMENSION AMAT(NDIM,NDIM) 1859C 1860 CALL QENTER('MATSYM') 1861 ISYM = 1 1862 IASYM = 2 1863 DO 220 J = 1,N 1864 DO 210 I = 1,J 1865 AMATS = ABS(AMAT(I,J) + AMAT(J,I)) 1866 AMATA = ABS(AMAT(I,J) - AMAT(J,I)) 1867 IF (AMATS .GT. THRZER) IASYM = 0 1868C ... i.e. not antisymmetric 1869 IF (AMATA .GT. THRZER) ISYM = 0 1870C ... i.e. not symmetric 1871 210 CONTINUE 1872 220 CONTINUE 1873 MATSYM = ISYM + IASYM 1874 CALL QEXIT('MATSYM') 1875 RETURN 1876 END 1877C /* Deck rewmot */ 1878 SUBROUTINE REWSPL(LU) 1879C 1880C Short interface routine for rewinding a file that may have been 1881C split, to ensure that we search for labels on the first of the split 1882C files. This routines preserves the UNIT number. 1883C 1884C K.Ruud, April 19 (2000) 1885C 1886#if !defined (VAR_SPLITFILES) 1887#include "priunit.h" 1888 INTEGER LU 1889 IF (LU .LT. 1) THEN 1890 WRITE (LUPRI,'(/A,I10)') 1891 & 'ERROR in REWSPL: negative unit number =',LU 1892 CALL QUIT('REWSPL called with negative unit number LU') 1893 END IF 1894 REWIND (LU) 1895#else 1896#include "implicit.h" 1897#include "dummy.h" 1898 CHARACTER FNNAME*80 1899C 1900 INQUIRE (UNIT=LU,NAME=FNNAME) 1901 LN = 1 1902 14 CONTINUE 1903 IF (FNNAME(LN:LN) .EQ. '-') THEN 1904 LN = LN - 1 1905 LUBKP = LU 1906 CALL GPCLOSE(LU,'KEEP') 1907 LU = LUBKP 1908 CALL GPOPEN(LU,FNNAME(1:LN),'OLD','SEQUENTIAL', 1909 & 'UNFORMATTED',IDUMMY,.FALSE.) 1910 INQUIRE (UNIT=LU,NAME=FNNAME) 1911 GOTO 15 1912 ELSE IF (FNNAME(LN:LN) .EQ. ' ') THEN 1913 GOTO 15 1914 END IF 1915 LN = LN + 1 1916 GOTO 14 1917 15 CONTINUE 1918 REWIND (LU) 1919#endif 1920 RETURN 1921 END 1922C /* Deck fndmin */ 1923 SUBROUTINE FNDMIN(NELMNT,IPLACE,VEC,NVEC,WRK,LWRK) 1924C 1925C 23-Nov-2000 hjaaj (FNDMIN = CIFNMN from 12-Aug-1990 hjaaj) 1926C (Find minimum elemnts) 1927C 1928C Purpose: 1929C Return in IPLACE addresses on lowest NELMNT elements in VEC. 1930C 1931#include "implicit.h" 1932 DIMENSION VEC(NVEC), IPLACE(NELMNT), WRK(LWRK) 1933C 1934 IF (LWRK .LT. NELMNT) THEN 1935 CALL QUIT('FNDMIN: Insufficient memory (LWRK .lt. NELMNT)') 1936 END IF 1937 IF (NELMNT .GT. NVEC) THEN 1938 CALL QUIT('FNDMIN ERROR: NELMNT .gt. NVEC') 1939 END IF 1940C 1941C Sort first NELMNT elements of VEC 1942C 1943 DO 120 I = 1,NELMNT 1944 VECI = VEC(I) 1945 DO 115 J = 1,(I-1) 1946 IF (WRK(J) .GT. VECI) THEN 1947 DO 112 K = I,(J+1),-1 1948 WRK(K) = WRK(K-1) 1949 IPLACE(K) = IPLACE(K-1) 1950 112 CONTINUE 1951 IPLACE(J) = I 1952 WRK(J) = VECI 1953 GO TO 120 1954 END IF 1955 115 CONTINUE 1956 IPLACE(I) = I 1957 WRK(I) = VECI 1958 120 CONTINUE 1959C 1960C Find lowest elements by insertion sort 1961C 1962 XHGH = WRK(NELMNT) 1963 DO 140 I = NELMNT+1,NVEC 1964 IF (VEC(I).GE.XHGH) GO TO 140 1965 DO 130 J = 1,NELMNT 1966 IF (VEC(I) .LT. WRK(J)) THEN 1967 DO 135 K = NELMNT,(J+1),-1 1968 WRK(K) = WRK(K-1) 1969 IPLACE(K) = IPLACE(K-1) 1970 135 CONTINUE 1971 IPLACE(J) = I 1972 WRK(J) = VEC(I) 1973 XHGH = WRK(NELMNT) 1974 GO TO 140 1975 END IF 1976 130 CONTINUE 1977 140 CONTINUE 1978 RETURN 1979 END 1980 SUBROUTINE COMPOFF(WRK, WORK, KBASE) 1981C 1982C Nov. 2004 Hans Joergen Aa. Jensen; revised Dec. 2010 1983C 1984C COMPOFF - Compute off-set from WORK to WRK 1985C i.e. WRK(I) = WORK(KBASE+I) 1986C 1987 REAL*8 :: WRK(*), WORK(*) 1988 integer(8), intent(inout) :: kbase 1989 integer(8) :: kaddr 1990 1991 !print *,'LOC(WRK) ',LOC(WRK) 1992 !print *,'LOC(WORK)',LOC(WORK) 1993 KADDR = LOC(WRK) - LOC(WORK) 1994 KBASE = KADDR / 8 + 1 ! from byte address to REAL*8 address 1995 1996 !print *,'COMPOFF: KADDR, KBASE', KADDR, KBASE, KBASE*8 1997 END 1998C -- end of gphjj.F -- 1999