1c $Id$ 2c 3c These are a bunch of utility routines that were in hessian/rhf_hessian.F 4c that were moved here since the hessian code no longer depends on them. 5c It may be that this needs to go to the util directory in the future, but 6c a lot of the functionality is not just utility routines. 7c 8 SUBROUTINE HND_REWFIL(NFT) 9 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 10 REWIND NFT 11 RETURN 12 END 13c 14 subroutine hnd_nwhnd_tran(dnwc,dhnd,ndim) 15 implicit double precision (a-h,o-z) 16#include "global.fh" 17 parameter (mxnbf =2048) 18 common/hnd_nwtohnd/inw_to_hnd(mxnbf) 19 common/hnd_iofile/ir,iw,ip 20 common/hnd_facntoh/fac_nwthnd(mxnbf) 21 dimension dnwc(ndim,ndim),dhnd(ndim,ndim) 22c 23 logical out 24 out =.false. 25c 26 if(out.and.ga_nodeid().eq.0) then 27 write(iw,*) ' in routine ... nwhnd_tran ... ' 28 endif 29c 30c ----- matrices tran. from nwchem to hondo shell order----- 31c 32 do j=1,ndim 33 do i=1,ndim 34 35 dhnd(inw_to_hnd(i),inw_to_hnd(j))= 36 & dnwc(i,j)/(fac_nwthnd(i)*fac_nwthnd(j)) 37 38 enddo 39 enddo 40c 41 return 42 end 43c 44 SUBROUTINE HND_DAREAD(IDAF,IODA,IX,NX,IDAR) 45 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 46 COMMON/HND_MACHIN/ISINGL,NBITS 47 COMMON/HND_DAFNAV/IXDDAF(2048),NAV10,NAV20 48 DIMENSION IODA(1),IX(1) 49 DIMENSION IXSDAF(1024) 50 EQUIVALENCE (IXSDAF(1),IXDDAF(1)) 51 DATA IXMAX /1024/ 52C 53 IF(IDAF.EQ.10) NAV10=IODA(IDAR+0) 54 IF(IDAF.EQ.20) NAV20=IODA(IDAR+0) 55 MAXIX=IXMAX*ISINGL 56 LDAR = NX*ISINGL 57C 58 MAX=0 59 10 MIN=MAX+1 60 MAX=MAX+MAXIX 61 IF(MAX.GT.LDAR) MAX=LDAR 62 IF(IDAF.EQ.10.AND.ISINGL.EQ.1) READ(IDAF,REC=NAV10) IXSDAF 63 IF(IDAF.EQ.10.AND.ISINGL.EQ.2) READ(IDAF,REC=NAV10) IXDDAF 64 IF(IDAF.EQ.20.AND.ISINGL.EQ.1) READ(IDAF,REC=NAV20) IXSDAF 65 IF(IDAF.EQ.20.AND.ISINGL.EQ.2) READ(IDAF,REC=NAV20) IXDDAF 66 DO 20 I=MIN,MAX 67 20 IX(I)=IXDDAF(I-MIN+1) 68 IF(IDAF.EQ.10) NAV10=NAV10+1 69 IF(IDAF.EQ.20) NAV20=NAV20+1 70 IF(MAX.LT.LDAR) GO TO 10 71 RETURN 72 END 73C 74 SUBROUTINE HND_DAWRIT(IDAF,IODA,IX,NX,IDAR,NAV) 75 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 76 CHARACTER*8 ERRMSG 77 PARAMETER (MXIODA=255) 78 LOGICAL OUT 79 COMMON/HND_IOFILE/IR,IW,IP 80 COMMON/HND_MACHIN/ISINGL,NBITS 81 COMMON/HND_DAFNAV/IXDDAF(2048),NAV10,NAV20 82 COMMON/HND_DAFREC/LENDAR(MXIODA) 83 DIMENSION ERRMSG(3) 84 DIMENSION IODA(1),IX(1) 85 DIMENSION IXSDAF(1024) 86 EQUIVALENCE (IXSDAF(1),IXDDAF(1)) 87 DATA ERRMSG /'PROGRAM ','STOP IN ','-DAWRIT-'/ 88 DATA IXMAX /1024/ 89C 90 OUT=.FALSE. 91C 92 MAXIX=IXMAX*ISINGL 93 LDAR = NX*ISINGL 94 IF(IODA(IDAR+0).NE.0) GO TO 100 95C 96C ----- FIRST WRITE ----- 97C 98 IODA(IDAR+0)=NAV 99 IF(IDAF.EQ.10) NAV10=NAV 100 IF(IDAF.EQ.20) NAV20=NAV 101C 102 IF(IDAF.EQ.10) THEN 103 LENDAR(IDAR+0)=NX 104 IF(OUT) THEN 105 WRITE(IW,9998) IDAR,NX 106 ENDIF 107 ENDIF 108C 109 MAX=0 110 10 MIN=MAX+1 111 MAX=MAX+MAXIX 112 IF(MAX.GT.LDAR) MAX=LDAR 113 DO 20 I=MIN,MAX 114 20 IXDDAF(I-MIN+1)=IX(I) 115 IF(IDAF.EQ.10.AND.ISINGL.EQ.1) WRITE(IDAF,REC=NAV10) IXSDAF 116 IF(IDAF.EQ.10.AND.ISINGL.EQ.2) WRITE(IDAF,REC=NAV10) IXDDAF 117 IF(IDAF.EQ.20.AND.ISINGL.EQ.1) WRITE(IDAF,REC=NAV20) IXSDAF 118 IF(IDAF.EQ.20.AND.ISINGL.EQ.2) WRITE(IDAF,REC=NAV20) IXDDAF 119 IF(IDAF.EQ.10) NAV10=NAV10+1 120 IF(IDAF.EQ.20) NAV20=NAV20+1 121 IF(MAX.LT.LDAR) GO TO 10 122C 123 IF(IDAF.EQ.10) NAV=NAV10 124 IF(IDAF.EQ.20) NAV=NAV20 125 RETURN 126C 127C ----- REWRITE ----- 128C 129 100 CONTINUE 130 IF(IDAF.EQ.10) NAV10=IODA(IDAR+0) 131 IF(IDAF.EQ.20) NAV20=IODA(IDAR+0) 132C 133 IF(IDAF.EQ.10) THEN 134 IF(OUT) THEN 135 WRITE(IW,9997) IDAR,NX 136 ENDIF 137 NX0=LENDAR(IDAR+0) 138 IF(NX.GT.NX0) THEN 139 IF(OUT) THEN 140 WRITE(IW,9999) IDAR,NX,NX0 141 CALL HND_HNDERR(3,ERRMSG) 142 ENDIF 143 ELSEIF(NX.LT.NX0) THEN 144 IF(OUT) THEN 145 WRITE(IW,9999) IDAR,NX,NX0 146 ENDIF 147 ENDIF 148 ENDIF 149C 150 MAX=0 151 110 MIN=MAX+1 152 MAX=MAX+MAXIX 153 IF(MAX.GT.LDAR) MAX=LDAR 154 DO 120 I=MIN,MAX 155 120 IXDDAF(I-MIN+1)=IX(I) 156 IF(IDAF.EQ.10.AND.ISINGL.EQ.1) WRITE(IDAF,REC=NAV10) IXSDAF 157 IF(IDAF.EQ.10.AND.ISINGL.EQ.2) WRITE(IDAF,REC=NAV10) IXDDAF 158 IF(IDAF.EQ.20.AND.ISINGL.EQ.1) WRITE(IDAF,REC=NAV20) IXSDAF 159 IF(IDAF.EQ.20.AND.ISINGL.EQ.2) WRITE(IDAF,REC=NAV20) IXDDAF 160 IF(IDAF.EQ.10) NAV10=NAV10+1 161 IF(IDAF.EQ.20) NAV20=NAV20+1 162 IF(MAX.LT.LDAR) GO TO 110 163C 164 RETURN 165 9999 FORMAT(' INCONSISTENT RECORD LENGTH FOR -IDAR- = ',I4,/, 166 1 ' NX,NX0 = ',2I10) 167 9998 FORMAT(' FIRST-WRITE, IDAR = ',I5,' LENGTH = ',I10) 168 9997 FORMAT(' RE-WRITE, IDAR = ',I5,' LENGTH = ',I10) 169 END 170C 171 subroutine hnd_hndnw_tran(dhnd,dnwc,ndim) 172 implicit double precision (a-h,o-z) 173#include "global.fh" 174 parameter (mxnbf =2048) 175 common/hnd_hndtonw/ihnd_to_nw(mxnbf) 176 common/hnd_iofile/ir,iw,ip 177 dimension dnwc(ndim,ndim),dhnd(ndim,ndim) 178 common/hnd_facntoh/fac_nwthnd(mxnbf) 179c 180 logical out 181 out =.false. 182c 183 if(out.and.ga_nodeid().eq.0) then 184 write(iw,*) ' in routine ... hndnw_tran ... ' 185 endif 186c 187c ----- matrices tran. from hondo to nwchem shell order----- 188c 189 do j=1,ndim 190 do i=1,ndim 191 dnwc(ihnd_to_nw(i),ihnd_to_nw(j))= 192 & dhnd(i,j)*(fac_nwthnd(ihnd_to_nw(i))* 193 & fac_nwthnd(ihnd_to_nw(j)) ) 194 enddo 195 enddo 196c 197 return 198 end 199c 200 SUBROUTINE HND_DSXYZ 201 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 202C 203C ----- GAUSS-HERMITE QUADRATURE USING MINIMUM POINT FORMULA ----- 204C 205#include "hnd_whermt.fh" 206c 207 COMMON/HND_XYZDER/XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ 208 1 ,NI,NJ 209 1 ,CX,CY,CZ 210 DIMENSION MIN(7),MAX(7) 211 DATA MIN /1,2,4, 7,11,16,22/ 212 DATA MAX /1,3,6,10,15,21,28/ 213 DATA ZERO /0.0D+00/ 214C 215 XINT=ZERO 216 YINT=ZERO 217 ZINT=ZERO 218 NPTS=(NI+NJ-2)/2+1 219 IMIN=MIN(NPTS) 220 IMAX=MAX(NPTS) 221 DO 16 I=IMIN,IMAX 222 DUM=W(I) 223 PX=DUM 224 PY=DUM 225 PZ=DUM 226 DUM=H(I)*T 227 PTX=DUM+X0 228 PTY=DUM+Y0 229 PTZ=DUM+Z0 230 AX=PTX-XI 231 AY=PTY-YI 232 AZ=PTZ-ZI 233 BX=PTX-XJ 234 BY=PTY-YJ 235 BZ=PTZ-ZJ 236 GO TO (7,6,5,4,3,2,1),NI 237 1 PX=PX*AX 238 PY=PY*AY 239 PZ=PZ*AZ 240 2 PX=PX*AX 241 PY=PY*AY 242 PZ=PZ*AZ 243 3 PX=PX*AX 244 PY=PY*AY 245 PZ=PZ*AZ 246 4 PX=PX*AX 247 PY=PY*AY 248 PZ=PZ*AZ 249 5 PX=PX*AX 250 PY=PY*AY 251 PZ=PZ*AZ 252 6 PX=PX*AX 253 PY=PY*AY 254 PZ=PZ*AZ 255 7 GO TO (15,14,13,12,11,10,9,8),NJ 256 8 PX=PX*BX 257 PY=PY*BY 258 PZ=PZ*BZ 259 9 PX=PX*BX 260 PY=PY*BY 261 PZ=PZ*BZ 262 10 PX=PX*BX 263 PY=PY*BY 264 PZ=PZ*BZ 265 11 PX=PX*BX 266 PY=PY*BY 267 PZ=PZ*BZ 268 12 PX=PX*BX 269 PY=PY*BY 270 PZ=PZ*BZ 271 13 PX=PX*BX 272 PY=PY*BY 273 PZ=PZ*BZ 274 14 PX=PX*BX 275 PY=PY*BY 276 PZ=PZ*BZ 277 15 CONTINUE 278 XINT=XINT+PX 279 YINT=YINT+PY 280 ZINT=ZINT+PZ 281 16 CONTINUE 282 RETURN 283 END 284C 285 SUBROUTINE HND_PREAD(IJK,XX,IX,NXX,MAX) 286 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 287 LOGICAL PACK2E 288 COMMON/HND_PCKLAB/LABSIZ 289 COMMON/HND_PCKOPT/NHEX,NTUPL,PACK2E 290 COMMON/HND_INTOPT/NOPK,NOK,NOSQUR 291 COMMON/HND_INTFIL/NINTMX 292 DIMENSION XX(MAX),IX(LABSIZ*MAX) 293 IF(PACK2E) GO TO 10 294 READ(IJK) XX,IX,NX 295 NXX=NX 296 RETURN 297 10 CALL HND_READPK(IJK,XX,XX, NXX,NH,IERR,IEND) 298 CALL HND_READPN(IJK,IX,IX,LABSIZ*NXX,NT,IERR,IEND) 299 RETURN 300 END 301C 302 SUBROUTINE HND_PRTRLAB(D,LAB,N) 303 IMPLICIT REAL*8 (A-H,O-Z) 304#include "nwc_const.fh" 305C 306C ----- PRINT OUT A TRIANGULAR MATRIX ----- 307C 308 COMMON/HND_IOFILE/IR,IW,IP 309 DIMENSION D(1),DD(10) 310 character*8 LAB 311 DIMENSION LAB(1) 312C 313 LIST=1 314 IF(LIST.EQ.1) MAX=7 315C 316 IMAX = 0 317 100 IMIN = IMAX+1 318 IMAX = IMAX+MAX 319 IF (IMAX .GT. N) IMAX = N 320 WRITE (IW,9008) 321 IF(LIST.EQ.1) WRITE (IW,9128) (I,I = IMIN,IMAX) 322 WRITE (IW,9008) 323 DO 160 J = 1,N 324 K = 0 325 DO 140 I = IMIN,IMAX 326 K = K+1 327 II = MAX0( I, J) 328 JJ = MIN0( I, J) 329 IJ = (II*(II-1))/2 + JJ 330 140 DD(K) = D(IJ) 331 IF(LIST.EQ.1) WRITE(IW,9148) J,LAB(J),(DD(I),I = 1,K) 332 160 CONTINUE 333 IF (IMAX .LT. N) GO TO 100 334 RETURN 335 9008 FORMAT(/) 336 9128 FORMAT(14X,7(6X,I3,6X)) 337 9148 FORMAT(I3,1X,A8,2X,7F15.10) 338 END 339C 340 SUBROUTINE HND_PRTRL(D,N) 341 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 342#include "nwc_const.fh" 343C 344C ----- PRINT OUT A TRIANGULAR MATRIX WITH LABELS 345C 346 PARAMETER (MXATOM=NW_MAX_ATOM) 347 PARAMETER (MXBFN=3072) 348 COMMON/HND_IOFILE/IR,IW,IP 349 COMMON/HND_LISTNG/LIST 350 DIMENSION D(1),DD(10) 351 IF(LIST.EQ.0) MAX=10 352 IF(LIST.EQ.1) MAX=7 353 IF(LIST.EQ.2) MAX=7 354 IMAX = 0 355 100 IMIN = IMAX+1 356 IMAX = IMAX+MAX 357 IF (IMAX .GT. N) IMAX = N 358 WRITE (IW,9008) 359c%% 360C IF(LIST.EQ.0) WRITE (IW,9028) (BFLAB(I),I = IMIN,IMAX) 361C IF(LIST.EQ.1) WRITE (IW,9128) (BFLAB(I),I = IMIN,IMAX) 362C IF(LIST.EQ.2) WRITE (IW,9228) (BFLAB(I),I = IMIN,IMAX) 363c%% 364 WRITE (IW,9008) 365 DO 160 J = 1,N 366 K = 0 367 DO 140 I = IMIN,IMAX 368 K = K+1 369 II = MAX0( I, J) 370 JJ = MIN0( I, J) 371 IJ = (II*(II-1))/2 + JJ 372 140 DD(K) = D(IJ) 373c%% 374C IF(LIST.EQ.0) WRITE (IW,9048) J,BFLAB(J),(DD(I),I = 1,K) 375C IF(LIST.EQ.1) WRITE (IW,9148) J,BFLAB(J),(DD(I),I = 1,K) 376C IF(LIST.EQ.2) WRITE (IW,9248) J,BFLAB(J),(DD(I),I = 1,K) 377c%% 378 160 CONTINUE 379 IF (IMAX .LT. N) GO TO 100 380 RETURN 381 9008 FORMAT(/) 382 9028 FORMAT(15X,10(2X,A8,1X)) 383 9048 FORMAT(I5,2X,A8,10F11.5) 384 9128 FORMAT(15X,7(4X,A8,3X)) 385 9148 FORMAT(I5,2X,A8,7F15.10) 386 9228 FORMAT(15X,7(4X,A8,3X)) 387 9248 FORMAT(I5,2X,A8,7E15.8) 388 END 389c 390 SUBROUTINE HND_SPRTR(D,N) 391 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 392C 393C ----- PRINT OUT A SQUARE MATRIX ----- 394C 395 COMMON/HND_IOFILE/IR,IW,IP 396 COMMON/HND_LISTNG/LIST 397 DIMENSION D(1),DD(10) 398C 399 LIST=1 400 IF(LIST.EQ.0) MAX=10 401 IF(LIST.EQ.1) MAX=7 402 IF(LIST.EQ.2) MAX=7 403C 404 IMAX = 0 405 100 IMIN = IMAX+1 406 IMAX = IMAX+MAX 407 IF (IMAX .GT. N) IMAX = N 408 WRITE (IW,9008) 409 IF(LIST.EQ.0) WRITE (IW,9028) (I,I = IMIN,IMAX) 410 IF(LIST.EQ.1) WRITE (IW,9128) (I,I = IMIN,IMAX) 411 IF(LIST.EQ.2) WRITE (IW,9228) (I,I = IMIN,IMAX) 412 WRITE (IW,9008) 413 DO 160 J = 1,N 414 K = 0 415 DO 140 I = IMIN,IMAX 416 K = K+1 417C II = MAX0( I, J) 418C JJ = MIN0( I, J) 419C IJ = (II*(II-1))/2 + JJ 420 IJ = (J-1)*N + I 421 140 DD(K) = D(IJ) 422 IF(LIST.EQ.0) WRITE (IW,9048) J,(DD(I),I = 1,K) 423 IF(LIST.EQ.1) WRITE (IW,9148) J,(DD(I),I = 1,K) 424 IF(LIST.EQ.2) WRITE (IW,9248) J,(DD(I),I = 1,K) 425 160 CONTINUE 426 IF (IMAX .LT. N) GO TO 100 427 RETURN 428 9008 FORMAT(/) 429 9028 FORMAT(6X,10(4X,I3,4X)) 430 9048 FORMAT(I5,1X,10F11.5) 431 9128 FORMAT(6X,7(6X,I3,6X)) 432 9148 FORMAT(I5,1X,7F15.10) 433 9228 FORMAT(6X,7(6X,I3,6X)) 434 9248 FORMAT(I5,1X,7E15.8) 435 END 436c 437 SUBROUTINE HND_SECOND(SEC,WSEC) 438 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 439 COMMON/HND_CLOCK0/SEC0,WSEC0 440 COMMON/HND_CLOCKS/DAT,TIM,CPUSEC,TOTSEC,DELCPU,DELTOT 441 DIMENSION DIAL(6) 442 EQUIVALENCE (DIAL(1),DAT) 443C 444C ----- RETURN ELAPSED - CPU- TIME IN SECONDS IN - SEC- ----- 445C ----- RETURN ELAPSED -WALL- TIME IN SECONDS IN -WSEC- ----- 446C 447 CALL HND_SYSCLK(DIAL) 448 SEC= CPUSEC 449 SEC= SEC- SEC0 450 WSEC= TOTSEC 451 WSEC=WSEC-WSEC0 452 RETURN 453 END 454C 455 SUBROUTINE HND_JKWRYS(RWV,ABV,NUMG) 456 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 457#include "hnd_rys.fh" 458 DIMENSION RWV(2,NUMG,1),ABV(5,1) 459C 460 IF(MROOTS.GT.5) GO TO 100 461C 462C ----- MROOTS .LE. 5 ----- 463C 464 DO 20 NG=1,NUMG 465 XX=ABV(5,NG) 466 IF(MROOTS.LE.3) CALL HND_RT123 467 IF(MROOTS.EQ.4) CALL HND_ROOT4 468 IF(MROOTS.EQ.5) CALL HND_ROOT5 469 DO 10 NR=1,MROOTS 470 RWV(1,NG,NR)=U(NR) 471 RWV(2,NG,NR)=W(NR) 472 10 CONTINUE 473 20 CONTINUE 474 RETURN 475C 476 100 CONTINUE 477C 478C ----- MROOTS .GT. 5 ----- 479C 480 DO 120 NG=1,NUMG 481 YY=ABV(5,NG) 482 CALL HND_DROOT 483 DO 110 NR=1,MROOTS 484 RWV(1,NG,NR)=U9(NR) 485 RWV(2,NG,NR)=W9(NR) 486 110 CONTINUE 487 120 CONTINUE 488 RETURN 489 END 490C 491 SUBROUTINE HND_JKBCDF(B00,B01,B10,C00,D00,F00,DIJ,DKL, 492 1 ABV,CV,RWV,NUMG,NROOTS) 493 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 494 LOGICAL NMAXS,NMAXP,MMAXS,MMAXP 495 COMMON/HND_SHLGNM/NMAXS,NMAXP,MMAXS,MMAXP 496 DIMENSION B00(NUMG,NROOTS,3),B01(NUMG,NROOTS,3),B10(NUMG,NROOTS,3) 497 DIMENSION C00(NUMG,NROOTS,3) 498 DIMENSION D00(NUMG,NROOTS,3) 499 DIMENSION F00(NUMG,NROOTS,3),DIJ(NUMG,NROOTS,3),DKL(NUMG,NROOTS,3) 500 DIMENSION ABV(5,1),CV(18,1) 501 DIMENSION RWV(2,NUMG,NROOTS) 502 DATA PT5,ONE /0.5D+00,1.0D+00/ 503C 504 DO 40 NR=1,NROOTS 505 DO 30 NG=1,NUMG 506 AA =ABV(1,NG) 507 BB =ABV(2,NG) 508 RHO=ABV(3,NG) 509 QAB=ABV(4,NG) 510 UU =RHO*RWV(1,NG,NR) 511 WW = RWV(2,NG,NR) 512 AAUU=AA+UU 513 BBUU=BB+UU 514 F00(NG,NR,1)=WW*QAB 515 F00(NG,NR,2)=ONE 516 F00(NG,NR,3)=ONE 517 DUM2=PT5/(AA*BB+UU*(AA+BB)) 518 AUDUM=AAUU*DUM2 519 BUDUM=BBUU*DUM2 520 UDUM= UU*DUM2 521 B00(NG,NR,1)= UDUM 522 B00(NG,NR,2)= UDUM 523 B00(NG,NR,3)= UDUM 524 B01(NG,NR,1)=AUDUM 525 B01(NG,NR,2)=AUDUM 526 B01(NG,NR,3)=AUDUM 527 B10(NG,NR,1)=BUDUM 528 B10(NG,NR,2)=BUDUM 529 B10(NG,NR,3)=BUDUM 530 UDUM= UDUM+ UDUM 531 IF(MMAXS) GO TO 10 532 AUDUM=AUDUM+AUDUM 533 D00(NG,NR,1)= UDUM*CV( 1,NG) + AUDUM*CV( 2,NG) 534 D00(NG,NR,2)= UDUM*CV( 3,NG) + AUDUM*CV( 4,NG) 535 D00(NG,NR,3)= UDUM*CV( 5,NG) + AUDUM*CV( 6,NG) 536 10 IF(NMAXS) GO TO 20 537 BUDUM=BUDUM+BUDUM 538 C00(NG,NR,1)= UDUM*CV( 8,NG) + BUDUM*CV( 7,NG) 539 C00(NG,NR,2)= UDUM*CV(10,NG) + BUDUM*CV( 9,NG) 540 C00(NG,NR,3)= UDUM*CV(12,NG) + BUDUM*CV(11,NG) 541 20 CONTINUE 542C 543 30 CONTINUE 544 40 CONTINUE 545C 546 DO 60 NR=1,NROOTS 547 DO 50 NG=1,NUMG 548 DIJ(NG,NR,1)=CV(13,NG) 549 DIJ(NG,NR,2)=CV(14,NG) 550 DIJ(NG,NR,3)=CV(15,NG) 551 DKL(NG,NR,1)=CV(16,NG) 552 DKL(NG,NR,2)=CV(17,NG) 553 DKL(NG,NR,3)=CV(18,NG) 554 50 CONTINUE 555 60 CONTINUE 556C 557 RETURN 558 END 559C 560 SUBROUTINE HND_JKGNMS(GNM,NG,NMAX,MMAX, 561 1 B00,B01,B10,C00,D00,F00) 562 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 563 LOGICAL NMAXS,NMAXP,MMAXS,MMAXP 564 COMMON/HND_SHLGNM/NMAXS,NMAXP,MMAXS,MMAXP 565 DIMENSION GNM(NG,NMAX,MMAX) 566 DIMENSION C00(NG),D00(NG),F00(NG) 567 DIMENSION B00(NG,1),B01(NG,1),B10(NG,1) 568C 569C ----- G(0,0) ----- 570C 571 DO 10 IG=1,NG 572 GNM(IG,1,1)=F00(IG) 573 10 CONTINUE 574 IF(NMAXS.AND.MMAXS) RETURN 575 IF(NMAXS) GO TO 30 576C 577C ----- G(1,0) = C00 * G(0,0) ----- 578C 579 DO 20 IG=1,NG 580 GNM(IG,2,1)=C00(IG)*GNM(IG,1,1) 581 20 CONTINUE 582C 583 30 CONTINUE 584 IF(MMAXS) GO TO 60 585C 586C ----- G(0,1) = D00 * G(0,0) ----- 587C 588 DO 40 IG=1,NG 589 GNM(IG,1,2)=D00(IG)*GNM(IG,1,1) 590 40 CONTINUE 591 IF(NMAXS) GO TO 60 592C 593C ----- G(1,1) = B00 * G(0,0) + D00 * G(1,0) ----- 594C 595 DO 50 IG=1,NG 596 GNM(IG,2,2)=B00(IG,1)*GNM(IG,1,1)+D00(IG)*GNM(IG,2,1) 597 50 CONTINUE 598C 599 60 CONTINUE 600 MAX=MAX0(NMAX-1,MMAX-1) 601 DO 70 M=2,MAX 602 DO 70 IG=1,NG 603 70 B00(IG,M)=B00(IG,M-1)+B00(IG,1) 604C 605 IF(NMAXP) GO TO 120 606C 607C ----- G(N+1,0) = N * B10 * G(N-1,0) + C00 * G(N,0) ----- 608C 609 DO 80 N=2,NMAX-1 610 DO 80 IG=1,NG 611 B10(IG,N)=B10(IG,N-1)+B10(IG,1) 612 80 CONTINUE 613 DO 90 N=2,NMAX-1 614 DO 90 IG=1,NG 615 GNM(IG,N+1,1)=B10(IG,N-1)*GNM(IG,N-1,1)+C00(IG)*GNM(IG,N,1) 616 90 CONTINUE 617 IF(MMAXS) GO TO 110 618C 619C ----- G(N,1) = N * B00 * G(N-1,0) + D00 * G(N,0) ----- 620C 621 DO 100 N=2,NMAX-1 622 DO 100 IG=1,NG 623 GNM(IG,N+1,2)=B00(IG,N)*GNM(IG,N,1)+D00(IG)*GNM(IG,N+1,1) 624 100 CONTINUE 625C 626 110 CONTINUE 627C 628 120 CONTINUE 629 IF(MMAXP) GO TO 170 630C 631C ----- G(0,M+1) = M * B01 * G(0,M-1) + D00 * G(O,M) ----- 632C 633 DO 130 M=2,MMAX-1 634 DO 130 IG=1,NG 635 B01(IG,M)=B01(IG,M-1)+B01(IG,1) 636 130 CONTINUE 637 DO 140 M=2,MMAX-1 638 DO 140 IG=1,NG 639 GNM(IG,1,M+1)=B01(IG,M-1)*GNM(IG,1,M-1)+D00(IG)*GNM(IG,1,M) 640 140 CONTINUE 641 IF(NMAXS) GO TO 160 642C 643C ----- G(1,M) = M * B00 * G(0,M-1) + C00 * G(0,M) ----- 644C 645 DO 150 M=2,MMAX-1 646 DO 150 IG=1,NG 647 GNM(IG,2,M+1)=B00(IG,M)*GNM(IG,1,M)+C00(IG)*GNM(IG,1,M+1) 648 150 CONTINUE 649C 650 160 CONTINUE 651 170 IF(NMAXP.OR.MMAXP) RETURN 652C 653C ----- G(N+1,M) = N * B10 * G(N-1,M ) ----- 654C + C00 * G(N ,M ) 655C + M * B00 * G(N ,M-1) 656C 657 DO 180 M=2,MMAX-1 658 DO 180 N=2,NMAX-1 659 DO 180 IG=1,NG 660 GNM(IG,N+1,M+1)=B10(IG,N-1)*GNM(IG,N-1,M+1)+ 661 1 C00(IG )*GNM(IG,N ,M+1)+ 662 2 B00(IG,M )*GNM(IG,N ,M ) 663 180 CONTINUE 664C 665 RETURN 666 END 667C 668 SUBROUTINE HND_JKGNMV(GNM,NG,NMAX,MMAX, 669 1 B00,B01,B10,C00,D00,F00) 670 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 671 LOGICAL NMAXS,NMAXP,MMAXS,MMAXP 672 COMMON/HND_SHLGNM/NMAXS,NMAXP,MMAXS,MMAXP 673 DIMENSION GNM(NG,NMAX,MMAX) 674 DIMENSION C00(NG),D00(NG),F00(NG) 675 DIMENSION B00(NG,1),B01(NG,1),B10(NG,1) 676C 677C ----- G(0,0) ----- 678C 679 DO 10 IG=1,NG 680 GNM(IG,1,1)=F00(IG) 681 10 CONTINUE 682 IF(NMAXS.AND.MMAXS) RETURN 683 IF(NMAXS) GO TO 30 684C 685C ----- G(1,0) = C00 * G(0,0) ----- 686C 687 DO 20 IG=1,NG 688 GNM(IG,2,1)=C00(IG)*GNM(IG,1,1) 689 20 CONTINUE 690C 691 30 CONTINUE 692 IF(MMAXS) GO TO 60 693C 694C ----- G(0,1) = D00 * G(0,0) ----- 695C 696 DO 40 IG=1,NG 697 GNM(IG,1,2)=D00(IG)*GNM(IG,1,1) 698 40 CONTINUE 699 IF(NMAXS) GO TO 60 700C 701C ----- G(1,1) = B00 * G(0,0) + D00 * G(1,0) ----- 702C 703 DO 50 IG=1,NG 704 GNM(IG,2,2)=B00(IG,1)*GNM(IG,1,1)+D00(IG)*GNM(IG,2,1) 705 50 CONTINUE 706C 707 60 CONTINUE 708 MAX=MAX0(NMAX-1,MMAX-1) 709 DO 70 IG=1,NG 710 DO 70 M=2,MAX 711 70 B00(IG,M)=B00(IG,M-1)+B00(IG,1) 712C 713 IF(NMAXP) GO TO 120 714C 715C ----- G(N+1,0) = N * B10 * G(N-1,0) + C00 * G(N,0) ----- 716C 717 DO 80 IG=1,NG 718 DO 80 N=2,NMAX-1 719 B10(IG,N)=B10(IG,N-1)+B10(IG,1) 720 80 CONTINUE 721 DO 90 IG=1,NG 722 DO 90 N=2,NMAX-1 723 GNM(IG,N+1,1)=B10(IG,N-1)*GNM(IG,N-1,1)+C00(IG)*GNM(IG,N,1) 724 90 CONTINUE 725 IF(MMAXS) GO TO 110 726C 727C ----- G(N,1) = N * B00 * G(N-1,0) + D00 * G(N,0) ----- 728C 729 DO 100 IG=1,NG 730 DO 100 N=2,NMAX-1 731 GNM(IG,N+1,2)=B00(IG,N)*GNM(IG,N,1)+D00(IG)*GNM(IG,N+1,1) 732 100 CONTINUE 733C 734 110 CONTINUE 735C 736 120 CONTINUE 737 IF(MMAXP) GO TO 170 738C 739C ----- G(0,M+1) = M * B01 * G(0,M-1) + D00 * G(O,M) ----- 740C 741 DO 130 IG=1,NG 742 DO 130 M=2,MMAX-1 743 B01(IG,M)=B01(IG,M-1)+B01(IG,1) 744 130 CONTINUE 745 DO 140 IG=1,NG 746 DO 140 M=2,MMAX-1 747 GNM(IG,1,M+1)=B01(IG,M-1)*GNM(IG,1,M-1)+D00(IG)*GNM(IG,1,M) 748 140 CONTINUE 749 IF(NMAXS) GO TO 160 750C 751C ----- G(1,M) = M * B00 * G(0,M-1) + C00 * G(0,M) ----- 752C 753 DO 150 IG=1,NG 754 DO 150 M=2,MMAX-1 755 GNM(IG,2,M+1)=B00(IG,M)*GNM(IG,1,M)+C00(IG)*GNM(IG,1,M+1) 756 150 CONTINUE 757C 758 160 CONTINUE 759 170 IF(NMAXP.OR.MMAXP) RETURN 760C 761C ----- G(N+1,M) = N * B10 * G(N-1,M ) ----- 762C + C00 * G(N ,M ) 763C + M * B00 * G(N ,M-1) 764C 765 DO 180 IG=1,NG 766 DO 180 N=2,NMAX-1 767 DO 180 M=2,MMAX-1 768 GNM(IG,N+1,M+1)=B10(IG,N-1)*GNM(IG,N-1,M+1)+ 769 1 C00(IG )*GNM(IG,N ,M+1)+ 770 2 B00(IG,M )*GNM(IG,N ,M ) 771 180 CONTINUE 772C 773 RETURN 774 END 775C 776 SUBROUTINE HND_JKXYZS(GIJKL,HIJKL,GNKL,HNKL,FNKL,GNM,HNM, 777 1 NG,NMAX,MMAX,NIMAX,NJMAX,NKMAX,NLMAX,DIJ,DKL, 778 2 EXPNDI,EXPNDK) 779 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 780 LOGICAL EXPNDI,EXPNDK 781 DIMENSION GIJKL(NG*NLMAX*NKMAX,NJMAX,NIMAX) 782 DIMENSION HIJKL(NG*NLMAX*NKMAX*NJMAX,NIMAX) 783 DIMENSION GNKL(NG,NLMAX,NKMAX,NMAX) 784 DIMENSION HNKL(NG*NLMAX*NKMAX,NMAX) 785 DIMENSION FNKL(NG*NLMAX*NKMAX*NMAX) 786 DIMENSION GNM(NG,NMAX,MMAX) 787 DIMENSION DIJ(NG) 788 DIMENSION DKL(NG) 789C 790C ----- G(N,K,L) ----- 791C 792 IF(EXPNDK) GO TO 40 793C 794 DO 30 NK=1,NKMAX 795 DO 10 NL=1,NLMAX 796 DO 10 N=1,NMAX 797 DO 10 IG=1,NG 798 10 GNKL(IG,NL,NK,N)=GNM(IG,N,NL) 799 IF(NK.EQ.NKMAX) GO TO 30 800 MAX=MMAX-NK 801 DO 20 M=1,MAX 802 DO 20 N=1,NMAX 803 DO 20 IG=1,NG 804 20 GNM(IG,N,M)=DKL(IG)*GNM(IG,N,M)+GNM(IG,N,M+1) 805 30 CONTINUE 806C 807 GO TO 100 808 40 CONTINUE 809C 810 DO 70 NL=1,NLMAX 811 DO 50 NK=1,NKMAX 812 DO 50 N=1,NMAX 813 DO 50 IG=1,NG 814 50 GNKL(IG,NL,NK,N)=GNM(IG,N,NK) 815 IF(NL.EQ.NLMAX) GO TO 70 816 MAX=MMAX-NL 817 DO 60 M=1,MAX 818 DO 60 N=1,NMAX 819 DO 60 IG=1,NG 820 60 GNM(IG,N,M)=DKL(IG)*GNM(IG,N,M)+GNM(IG,N,M+1) 821 70 CONTINUE 822C 823 100 CONTINUE 824C 825C ----- G(I,J,K,L) ----- 826C 827 IF(EXPNDI) GO TO 140 828C 829 DO 130 NI=1,NIMAX 830 DO 110 IGLKJ=1,NG*NLMAX*NKMAX*NJMAX 831 110 HIJKL(IGLKJ,NI)=FNKL(IGLKJ) 832 IF(NI.EQ.NIMAX) GO TO 130 833 MAX=NMAX-NI 834 DO 120 N=1,MAX 835 DO 120 NK=1,NKMAX 836 DO 120 NL=1,NLMAX 837 DO 120 IG=1,NG 838 120 GNKL(IG,NL,NK,N)=DIJ(IG)*GNKL(IG,NL,NK,N)+GNKL(IG,NL,NK,N+1) 839 130 CONTINUE 840C 841 RETURN 842 140 CONTINUE 843C 844 DO 170 NJ=1,NJMAX 845 DO 150 NI=1,NIMAX 846 DO 150 IGLK=1,NG*NLMAX*NKMAX 847 150 GIJKL(IGLK,NJ,NI)=HNKL(IGLK,NI) 848 IF(NJ.EQ.NJMAX) GO TO 170 849 MAX=NMAX-NJ 850 DO 160 N=1,MAX 851 DO 160 NK=1,NKMAX 852 DO 160 NL=1,NLMAX 853 DO 160 IG=1,NG 854 160 GNKL(IG,NL,NK,N)=DIJ(IG)*GNKL(IG,NL,NK,N)+GNKL(IG,NL,NK,N+1) 855 170 CONTINUE 856C 857 RETURN 858 END 859C 860 SUBROUTINE HND_JKXYZV(GIJKL,HIJKL,GNKL,HNKL,FNKL,GNM,HNM, 861 1 NG,NMAX,MMAX,NIMAX,NJMAX,NKMAX,NLMAX,DIJ,DKL, 862 2 EXPNDI,EXPNDK) 863 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 864 LOGICAL EXPNDI,EXPNDK 865 DIMENSION GIJKL(NG*NLMAX*NKMAX,NJMAX,NIMAX) 866 DIMENSION HIJKL(NG*NLMAX*NKMAX*NJMAX,NIMAX) 867 DIMENSION GNKL(NG,NLMAX,NKMAX,NMAX) 868 DIMENSION HNKL(NG*NLMAX*NKMAX,NMAX) 869 DIMENSION FNKL(NG*NLMAX*NKMAX*NMAX) 870 DIMENSION GNM(NG,NMAX,MMAX) 871 DIMENSION DIJ(NG) 872 DIMENSION DKL(NG) 873C 874C ----- G(N,K,L) ----- 875C 876 IF(EXPNDK) GO TO 40 877C 878 DO 30 NK=1,NKMAX 879 DO 10 IG=1,NG 880 DO 10 NL=1,NLMAX 881 DO 10 N=1,NMAX 882 10 GNKL(IG,NL,NK,N)=GNM(IG,N,NL) 883 IF(NK.EQ.NKMAX) GO TO 30 884 MAX=MMAX-NK 885 DO 20 IG=1,NG 886 DO 20 M=1,MAX 887 DO 20 N=1,NMAX 888 20 GNM(IG,N,M)=DKL(IG)*GNM(IG,N,M)+GNM(IG,N,M+1) 889 30 CONTINUE 890C 891 GO TO 100 892 40 CONTINUE 893C 894 DO 70 NL=1,NLMAX 895 DO 50 IG=1,NG 896 DO 50 NK=1,NKMAX 897 DO 50 N=1,NMAX 898 50 GNKL(IG,NL,NK,N)=GNM(IG,N,NK) 899 IF(NL.EQ.NLMAX) GO TO 70 900 MAX=MMAX-NL 901 DO 60 IG=1,NG 902 DO 60 N=1,NMAX 903 DO 60 M=1,MAX 904 60 GNM(IG,N,M)=DKL(IG)*GNM(IG,N,M)+GNM(IG,N,M+1) 905 70 CONTINUE 906C 907 100 CONTINUE 908C 909C ----- G(I,J,K,L) ----- 910C 911 IF(EXPNDI) GO TO 140 912C 913 DO 130 NI=1,NIMAX 914 DO 110 IGLKJ=1,NG*NLMAX*NKMAX*NJMAX 915 110 HIJKL(IGLKJ,NI)=FNKL(IGLKJ) 916 IF(NI.EQ.NIMAX) GO TO 130 917 MAX=NMAX-NI 918 DO 120 IG=1,NG 919 DO 120 NL=1,NLMAX 920 DO 120 NK=1,NKMAX 921 DO 120 N=1,MAX 922 120 GNKL(IG,NL,NK,N)=DIJ(IG)*GNKL(IG,NL,NK,N)+GNKL(IG,NL,NK,N+1) 923 130 CONTINUE 924C 925 RETURN 926 140 CONTINUE 927C 928 DO 170 NJ=1,NJMAX 929 DO 150 IGLK=1,NG*NLMAX*NKMAX 930 DO 150 NI=1,NIMAX 931 150 GIJKL(IGLK,NJ,NI)=HNKL(IGLK,NI) 932 IF(NJ.EQ.NJMAX) GO TO 170 933 MAX=NMAX-NJ 934 DO 160 IG=1,NG 935 DO 160 NL=1,NLMAX 936 DO 160 NK=1,NKMAX 937 DO 160 N=1,MAX 938 160 GNKL(IG,NL,NK,N)=DIJ(IG)*GNKL(IG,NL,NK,N)+GNKL(IG,NL,NK,N+1) 939 170 CONTINUE 940C 941 RETURN 942 END 943C 944 SUBROUTINE HND_READPK(IS,XX,YY,NXX,NH,IERR,IEND) 945 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 946 CHARACTER*8 ERRMSG 947 DIMENSION ERRMSG(3) 948 DATA ERRMSG /'PROGRAM ','STOP IN ','-READPK-'/ 949 CALL HND_HNDERR(3,ERRMSG) 950 RETURN 951 END 952c 953 SUBROUTINE HND_READPN(IS,IX,IY,NNX,NT,IERR,IEND) 954 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 955 CHARACTER*8 ERRMSG 956 DIMENSION ERRMSG(3) 957 DATA ERRMSG /'PROGRAM ','STOP IN ','-READPN-'/ 958 CALL HND_HNDERR(3,ERRMSG) 959 RETURN 960 END 961c 962 SUBROUTINE HND_SYSCLK(DIAL) 963 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 964 PARAMETER (TEN3=1.0D+03) 965C 966C ----- -CPU- AND -WALL- TIME ----- 967C 968C DIAL(1) = DATE 969C DIAL(2) = TIME 970C DIAL(3) = CPUSEC 971C DIAL(4) = TOTSEC 972C DIAL(5) = DELCPU 973C DIAL(6) = DELTOT 974C 975 DIMENSION DIAL(6) 976* REAL*4 CPUTIM,ETIME_ 977* TYPE TB_TYPE 978* SEQUENCE 979* REAL*4 USRTIM 980* REAL*4 SYSTIM 981* END TYPE 982c TYPE (TB_TYPE) ETIME_STRUCT 983c CPUTIM =ETIME_(ETIME_STRUCT) 984c DIAL(3)=CPUTIM 985* TOTTIM =TIMEF()/TEN3 986* DIAL(4)=TOTTIM 987 RETURN 988 END 989