1C FILE : printpkg.F 2C 3C... Dalton, pdpack/printpkg.F 4C... 5C... These routines are in the public domain and can be 6C... used freely in other programs. 7C... 8C 9C FILE: printpkg.F 10C 11C ==================================================================================== 12C For all OUT* routines, the NCTL parameter is used to select different output formats. 13C ------------------------------------------------------------------------------------ 14C NCTL > 0: fit in 80 columns 15C NCTL < 0: fit in 132 columns 16C abs(NCTL) .ge. 10: also print zero rows 17C abs(NCTL) =2 or =3: write old obsolete line printer carriage control characters 18C for double or triple space between lines 19C ==================================================================================== 20C 21C /* Deck outpak */ 22 SUBROUTINE OUTPAK (AMATRX,NROW,NCTL,LUPRI) 23C....................................................................... 24C Revised 16-Dec-1983 by Hans Jorgen Aa. Jensen. 25C 16-Jun-1986 hjaaj ( removed Hollerith ) 26C 27C OUTPAK PRINTS A REAL SYMMETRIC MATRIX STORED IN ROW-PACKED LOWER 28C 29C TRIANGULAR FORM (SEE DIAGRAM BELOW) IN FORMATTED FORM WITH NUMBERED 30C 31C ROWS AND COLUMNS. THE INPUT IS AS FOLLOWS: 32C 33C AMATRX(*)...........PACKED MATRIX 34C 35C NROW................NUMBER OF ROWS TO BE OUTPUT 36C 37C NCTL................CARRIAGE CONTROL FLAG: 1 FOR SINGLE SPACE, 38C 2 FOR DOUBLE SPACE, 39C 3 FOR TRIPLE SPACE. 40C 41C THE MATRIX ELEMENTS ARE ARRANGED IN STORAGE AS FOLLOWS: 42C 43C 1 44C 2 3 45C 4 5 6 46C 7 8 9 10 47C 11 12 13 14 15 48C 16 17 18 19 20 21 49C 22 23 24 25 26 27 28 50C AND SO ON. 51C 52C OUTPAK IS SET UP TO HANDLE 6 COLUMNS/PAGE WITH A 6F20.14 FORMAT 53C FOR THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, CHANGE 54C FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER OF 55C COLUMNS. 56C 57C AUTHOR: NELSON H.F. BEEBE, QUANTUM THEORY PROJECT, UNIVERSITY OF 58C FLORIDA, GAINESVILLE 59C..........VERSION = 09/05/73/03 60C....................................................................... 61C 62#include "implicit.h" 63 DIMENSION AMATRX(*) 64 INTEGER BEGIN 65 CHARACTER*1 ASA(3),BLANK,CTL 66 CHARACTER PFMT*20, COLUMN*8 67 PARAMETER (ZERO=0.D00, KCOLP=4, KCOLN=6) 68 PARAMETER (FFMIN=1.D-3, FFMAX = 1.D3) 69 LOGICAL IS_NAN 70 DATA COLUMN/'Column '/, ASA/' ', '0', '-'/, BLANK/' '/ 71C 72 IF (NCTL .LT. 0) THEN 73 KCOL = KCOLN 74 ELSE 75 KCOL = KCOLP 76 END IF 77 MCTL = ABS(NCTL) 78 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 79 CTL = ASA(MCTL) 80 ELSE 81 CTL = BLANK 82 END IF 83C 84 J = NROW*(NROW+1)/2 85 AMAX = ZERO 86 BMAX = ZERO 87 CMAX = ZERO 88 N_NAN = 0 89 DO I = 1,J 90 IF (IS_NAN(AMATRX(I),AMATRX(I))) THEN 91 N_NAN = N_NAN + 1 92 ELSE 93 AMAX = MAX( AMAX, ABS(AMATRX(I)) ) 94 BMAX = MAX( BMAX, AMATRX(I)) 95 CMAX = MIN( CMAX, AMATRX(I)) 96 END IF 97 END DO 98 IF (N_NAN .GT. 0) WRITE (LUPRI,'(/T6,A,I10,A)') 99 & 'WARNING: matrix contains',N_NAN,' NaN.' 100 IF (AMAX .EQ. ZERO) THEN 101 WRITE (LUPRI,'(/T6,A)') 'Zero matrix.' 102 GO TO 200 103 END IF 104 IF (FFMIN .LE. AMAX .AND. AMAX .LE. FFMAX) THEN 105C use F output format 106 PFMT = '(A1,I7,2X,8F15.8)' 107 ELSE 108C use 1PD output format 109 PFMT = '(A1,I7,2X,1P,8E15.6)' 110 END IF 111C 112C LAST IS THE LAST COLUMN NUMBER IN THE ROW CURRENTLY BEING PRINTED 113C 114 LAST = MIN(NROW,KCOL) 115C 116C BEGIN IS THE FIRST COLUMN NUMBER IN THE ROW CURRENTLY BEING PRINTED. 117C 118C.....BEGIN NON STANDARD DO LOOP. 119 BEGIN= 1 120 1050 NCOL = 1 121 WRITE (LUPRI,1000) (COLUMN,I,I = BEGIN,LAST) 122 DO 40 K = BEGIN,NROW 123 KTOTAL = (K*(K-1))/2 + BEGIN - 1 124 IF (MCTL .GE. 10) GO TO 20 ! also print zero rows 125 DO 10 I = 1,NCOL 126 IF (AMATRX(KTOTAL+I) .NE. ZERO) GO TO 20 127 10 CONTINUE 128 GO TO 30 129 20 WRITE (LUPRI,PFMT) CTL,K,(AMATRX(J+KTOTAL),J=1,NCOL) 130 30 IF (K .LT. (BEGIN+KCOL-1)) NCOL = NCOL + 1 131 40 CONTINUE 132 LAST = MIN(LAST+KCOL,NROW) 133 BEGIN= BEGIN + NCOL 134 IF (BEGIN.LE.NROW) GO TO 1050 135 200 CONTINUE 136 WRITE(LUPRI,'(A)') ' ==== End of matrix output ====' 137 RETURN 138C 139 1000 FORMAT (/10X,8(5X,A6,I4)) 140 END 141C /* Deck outpkb */ 142 SUBROUTINE OUTPKB (AMATRX,NDIM,NBLK,NCTL,LUPRI) 143C....................................................................... 144C 145C OUTPKB is OUTPAK modified for blocking as specified by 146C NDIM(NBLK). 16-Dec-1983 Hans Jorgen Aa. Jensen. 147C 148C 16-Jun-1986 hjaaj ( removed Hollerith ) 149C 150C OUTPAK PRINTS A REAL OMSYMMETRIC MATRIX STORED IN ROW-PACKED LOWER 151C TRIANGULAR FORM (SEE DIAGRAM BELOW) IN FORMATTED FORM WITH NUMBERED 152C ROWS AND COLUMNS. THE INPUT IS AS FOLLOWS: 153C 154C AMATRX(')...........PACKED MATRIX, blocked 155C 156C NDIM(').............dimension of each block 157C 158C NBLK................number of blocks 159C 160C NCTL................CARRIAGE CONTROL FLAG: 1 FOR SINGLE SPACE, 161C 2 FOR DOUBLE SPACE, 162C 3 FOR TRIPLE SPACE. 163C 164C THE MATRIX ELEMENTS in a block ARE ARRANGED IN STORAGE AS 165C FOLLOWS: 166C 167C 1 168C 2 3 169C 4 5 6 170C 7 8 9 10 171C 11 12 13 14 15 172C 16 17 18 19 20 21 173C 22 23 24 25 26 27 28 174C AND SO ON. 175C 176C OUTPAK IS SET UP TO HANDLE 6 COLUMNS/PAGE WITH A 6F20.14 FORMAT 177C FOR THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, CHANGE 178C FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER OF 179C COLUMNS. 180C 181C AUTHOR: NELSON H.F. BEEBE, QUANTUM THEORY PROJECT, UNIVERSITY OF 182C FLORIDA, GAINESVILLE 183C..........OUTPAK VERSION = 09/05/73/03 184C....................................................................... 185C 186#include "implicit.h" 187 INTEGER NDIM(NBLK),BEGIN 188 DIMENSION AMATRX(*) 189 CHARACTER*1 ASA(3),BLANK,CTL 190 CHARACTER PFMT*20, COLUMN*8 191 LOGICAL IS_NAN 192 PARAMETER (ZERO = 0.0D0, KCOLP=4, KCOLN=6) 193 PARAMETER (FFMIN=1.D-3, FFMAX = 1.D3) 194 DATA COLUMN/'Column '/, ASA/' ', '0', '-'/, BLANK/' '/ 195C 196 IF (NCTL .LT. 0) THEN 197 KCOL = KCOLN 198 ELSE 199 KCOL = KCOLP 200 END IF 201 MCTL = ABS(NCTL) 202 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 203 CTL = ASA(MCTL) 204 ELSE 205 CTL = BLANK 206 END IF 207C 208 MATLN = 0 209 DO 200 IBLK = 1,NBLK 210 MATLN = MATLN + NDIM(IBLK)*(NDIM(IBLK)+1)/2 211 200 CONTINUE 212C 213 AMAX = ZERO 214 N_NAN = 0 215 DO I=1,MATLN 216 IF ( IS_NAN(AMATRX(I),AMATRX(I)) ) THEN 217 N_NAN = N_NAN + 1 218 ELSE 219 AMAX = MAX( AMAX, ABS(AMATRX(I)) ) 220 END IF 221 END DO 222 IF (N_NAN .GT. 0) WRITE (LUPRI,'(/T6,A,I10,A)') 223 & 'WARNING: matrix contains',N_NAN,' NaN.' 224 IF (AMAX .EQ. ZERO) THEN 225 WRITE (LUPRI,3000) NBLK 226 GO TO 800 227 END IF 228 IF (FFMIN .LE. AMAX .AND. AMAX .LE. FFMAX) THEN 229C use F output format 230 PFMT = '(A1,I7,2X,8F15.8)' 231 ELSE 232C use 1PD output format 233 PFMT = '(A1,I7,2X,1P,8D15.6)' 234 END IF 235C 236 IOFF = 0 237 DO 100 IBLK = 1,NBLK 238 IDIM = NDIM(IBLK) 239 IF (IDIM.EQ.0) GO TO 100 240 IIDIM = IDIM*(IDIM+1)/2 241C 242 DO 5 I=1,IIDIM 243 IF (AMATRX(IOFF+I).NE.ZERO) GO TO 15 244 5 CONTINUE 245 WRITE (LUPRI,1100) IBLK 246 GO TO 90 247 15 CONTINUE 248 WRITE (LUPRI,1200) IBLK 249C 250C LAST IS THE LAST COLUMN NUMBER IN THE ROW CURRENTLY BEING PRINTED 251C 252 LAST = MIN(IDIM,KCOL) 253C 254C BEGIN IS THE FIRST COLUMN NUMBER IN THE ROW CURRENTLY BEING PRINTED. 255C 256 BEGIN = 1 257C.....BEGIN NON STANDARD DO LOOP. 258 1050 NCOL = 1 259 WRITE (LUPRI,1000) (COLUMN,I,I = BEGIN,LAST) 260 KTOTAL = IOFF + BEGIN*(BEGIN+1)/2 - 1 261 DO 40 K = BEGIN,IDIM 262 IF (MCTL .GE. 10) GO TO 20 ! also print zero rows 263 DO 10 I = 1,NCOL 264 IF (AMATRX(KTOTAL+I) .NE. ZERO) GO TO 20 265 10 CONTINUE 266 GO TO 30 267 20 WRITE (LUPRI,PFMT) CTL,K,(AMATRX(KTOTAL+J),J=1,NCOL) 268 30 IF (K .LT. (BEGIN+KCOL-1)) NCOL = NCOL + 1 269 KTOTAL = KTOTAL + K 270 40 CONTINUE 271 LAST = MIN(LAST+KCOL,IDIM) 272 BEGIN=BEGIN+NCOL 273 IF (BEGIN.LE.IDIM) GO TO 1050 274C 275 90 IOFF = IOFF + IIDIM 276 100 CONTINUE 277C 278 800 CONTINUE 279 WRITE(LUPRI,'(A)') ' ==== End of matrix output ====' 280 RETURN 281 3000 FORMAT (/5X,'All',I3,' blocks zero matrices.') 282 1100 FORMAT (/5X,'*** Block',I3,' zero matrix. ***') 283 1200 FORMAT (/5X,'*** Block',I3,' ***') 284 1000 FORMAT (/10X,8(5X,A6,I4)) 285 END 286C /* Deck outpks */ 287 SUBROUTINE OUTPKS (AMATRX,NDIM,NBLK,IREPO,NCTL,LUPRI) 288C....................................................................... 289C 290C OUTPKS is OUTPAK modified for blocking as specified by 291C NDIM(NBLK). 16-Dec-1983 Hans Jorgen Aa. Jensen. 292C 14-Dec-1988 tuh , accepts not totally symmetric elements 293C 294C 16-Jun-1986 hjaaj ( removed Hollerith ) 295C 296C OUTPAK PRINTS A REAL OMSYMMETRIC MATRIX STORED IN ROW-PACKED LOWER 297C TRIANGULAR FORM (SEE DIAGRAM BELOW) IN FORMATTED FORM WITH NUMBERED 298C ROWS AND COLUMNS. THE INPUT IS AS FOLLOWS: 299C 300C AMATRX(')...........PACKED MATRIX, blocked 301C 302C NDIM(').............dimension of each block 303C 304C NBLK................number of blocks 305C 306C NCTL................CARRIAGE CONTROL FLAG: 1 FOR SINGLE SPACE, 307C 2 FOR DOUBLE SPACE, 308C 3 FOR TRIPLE SPACE. 309C 310C THE MATRIX ELEMENTS in a block ARE ARRANGED IN STORAGE AS 311C FOLLOWS: 312C 313C 1 314C 2 3 315C 4 5 6 316C 7 8 9 10 317C 11 12 13 14 15 318C 16 17 18 19 20 21 319C 22 23 24 25 26 27 28 320C AND SO ON. 321C 322C OUTPAK IS SET UP TO HANDLE 6 COLUMNS/PAGE WITH A 6F20.14 FORMAT 323C FOR THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, CHANGE 324C FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER OF 325C COLUMNS. 326C 327C AUTHOR: NELSON H.F. BEEBE, QUANTUM THEORY PROJECT, UNIVERSITY OF 328C FLORIDA, GAINESVILLE 329C..........OUTPAK VERSION = 09/05/73/03 330C....................................................................... 331C 332#include "implicit.h" 333 INTEGER NDIM(NBLK),BEGIN 334 DIMENSION AMATRX(*) 335 CHARACTER*1 ASA(3),BLANK,CTL 336 CHARACTER PFMT*20, COLUMN*8 337 PARAMETER (D0 = 0.0D0, KCOLP=4, KCOLN=6) 338 PARAMETER (FFMIN=1.D-3, FFMAX = 1.D3) 339 DATA COLUMN/'Column '/, ASA/' ', '0', '-'/, BLANK/' '/ 340 341C 342 IF (IREPO .NE. 1) THEN 343 IBLK = 1 344 DO 500 IREPB = 0, NBLK - 1 345 IREPA = IEOR(IREPO - 1,IREPB) 346 IF (IREPA .GT. IREPB) THEN 347 NDIMA = NDIM(IREPA + 1) 348 NDIMB = NDIM(IREPB + 1) 349 WRITE (LUPRI,'(/A,2I5)') ' Symmetries: ',IREPA+1,IREPB+1 350 WRITE (LUPRI,'(/A,2I5)') ' Dimensions: ',NDIMA,NDIMB 351 IF (NDIMA.GT.0 .AND. NDIMB.GT.0) THEN 352 CALL OUTPUT(AMATRX(IBLK),1,NDIMA,1,NDIMB, 353 * NDIMA,NDIMB,NCTL,LUPRI) 354 IBLK = IBLK + NDIMA*NDIMB 355 END IF 356 ENDIF 357 500 CONTINUE 358 GO TO 800 359 END IF 360 IF (NCTL .LT. 0) THEN 361 KCOL = KCOLN 362 ELSE 363 KCOL = KCOLP 364 END IF 365 MCTL = ABS(NCTL) 366 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 367 CTL = ASA(MCTL) 368 ELSE 369 CTL = BLANK 370 END IF 371C 372 MATLN = 0 373 DO 200 IBLK = 1,NBLK 374 MATLN = MATLN + NDIM(IBLK)*(NDIM(IBLK)+1)/2 375 200 CONTINUE 376C 377 AMAX = D0 378 DO 205 I=1,MATLN 379 AMAX = MAX( AMAX, ABS(AMATRX(I)) ) 380 205 CONTINUE 381 IF (AMAX .EQ. D0) THEN 382 WRITE (LUPRI,3000) NBLK 383 GO TO 800 384 END IF 385 IF (FFMIN .LE. AMAX .AND. AMAX .LE. FFMAX) THEN 386C use F output format 387 PFMT = '(A1,I7,2X,8F15.8)' 388 ELSE 389C use 1PD output format 390 PFMT = '(A1,I7,2X,1P,8D15.6)' 391 END IF 392C 393 IOFF = 0 394 DO 100 IBLK = 1,NBLK 395 IDIM = NDIM(IBLK) 396 IF (IDIM.EQ.0) GO TO 100 397 IIDIM = IDIM*(IDIM+1)/2 398C 399 DO 5 I=1,IIDIM 400 IF (AMATRX(IOFF+I).NE.D0) GO TO 15 401 5 CONTINUE 402 WRITE (LUPRI,1100) IBLK 403 GO TO 90 404 15 CONTINUE 405 WRITE (LUPRI,1200) IBLK 406C 407C LAST IS THE LAST COLUMN NUMBER IN THE ROW CURRENTLY BEING PRINTED 408C 409 LAST = MIN(IDIM,KCOL) 410C 411C BEGIN IS THE FIRST COLUMN NUMBER IN THE ROW CURRENTLY BEING PRINTED. 412C 413 BEGIN = 1 414C.....BEGIN NON STANDARD DO LOOP. 415 1050 NCOL = 1 416 WRITE (LUPRI,1000) (COLUMN,I,I = BEGIN,LAST) 417 KTOTAL = IOFF + BEGIN*(BEGIN+1)/2 - 1 418 DO 40 K = BEGIN,IDIM 419 IF (MCTL .GE. 10) GO TO 20 ! also print zero rows 420 DO 10 I = 1,NCOL 421 IF (AMATRX(KTOTAL+I) .NE. D0) GO TO 20 422 10 CONTINUE 423 GO TO 30 424 20 WRITE (LUPRI,PFMT) CTL,K,(AMATRX(KTOTAL+J),J=1,NCOL) 425 30 IF (K .LT. (BEGIN+KCOL-1)) NCOL = NCOL + 1 426 KTOTAL = KTOTAL + K 427 40 CONTINUE 428 LAST = MIN(LAST+KCOL,IDIM) 429 BEGIN=BEGIN+NCOL 430 IF (BEGIN.LE.IDIM) GO TO 1050 431C 432 90 IOFF = IOFF + IIDIM 433 100 CONTINUE 434C 435 800 CONTINUE 436 WRITE(LUPRI,'(A)') ' ==== End of matrix output ====' 437 RETURN 438 3000 FORMAT (/5X,'All',I3,' blocks zero matrices.') 439 1100 FORMAT (/5X,'*** Block',I3,' zero matrix. ***') 440 1200 FORMAT (/5X,'*** Block',I3,' ***') 441 1000 FORMAT (/10X,8(5X,A6,I4)) 442 END 443C /* Deck outptb */ 444 SUBROUTINE OUTPTB (AMATRX,NDIM,NBLK,NCTL,LUPRI) 445C....................................................................... 446C 447C OUTPTB is OUTPUT modified for blocking as specified by 448C NDIM(NBLK). 16-Dec-1983 Hans Jorgen Aa. Jensen. 449C 450C 16-Jun-1986 hjaaj ( removed Hollerith ) 451C 452C OUTPUT PRINTS A REAL MATRIX IN FORMATTED FORM WITH NUMBERED ROWS 453C AND COLUMNS. THE INPUT IS AS FOLLOWS; 454C 455C AMATRX(',').........MATRIX TO BE OUTPUT, blocked 456C 457C NDIM(').............dimension of each block 458C 459C NBLK................number of blocks 460C 461C NCTL................CARRIAGE CONTROL FLAG; 1 FOR SINGLE SPACE 462C 2 FOR DOUBLE SPACE 463C 3 FOR TRIPLE SPACE 464C 465C THE PARAMETERS THAT FOLLOW MATRIX ARE ALL OF TYPE INTEGER. THE 466C PROGRAM IS SET UP TO HANDLE 5 COLUMNS/PAGE WITH A 1P,5D24.15 FORMAT 467C FOR 1THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, 468C CHANGE FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER 469C OF COLUMNS. 470C 471C AUTHOR; NELSON H.F. BEEBE, QUANTUM THEORY PROJECT, UNIVERSITY OF 472C FLORIDA, GAINESVILLE 473C REVISED; FEBRUARY 26, 1971 474C 475C....................................................................... 476C 477#include "implicit.h" 478 DIMENSION AMATRX(*) 479 INTEGER NDIM(NBLK), BEGIN, KCOL 480 CHARACTER*1 ASA(3), BLANK, CTL 481 CHARACTER PFMT*20, COLUMN*8 482 LOGICAL IS_NAN 483 PARAMETER (ZERO=0.D00, KCOLP=4, KCOLN=6) 484 PARAMETER (FFMIN=1.D-3, FFMAX = 1.D3) 485 DATA COLUMN/'Column '/, ASA/' ', '0', '-'/, BLANK/' '/ 486C 487 IF (NCTL .LT. 0) THEN 488 KCOL = KCOLN 489 ELSE 490 KCOL = KCOLP 491 END IF 492 MCTL = ABS(NCTL) 493 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 494 CTL = ASA(MCTL) 495 ELSE 496 CTL = BLANK 497 END IF 498C 499 MATDIM = 0 500 DO 200 IBLK = 1,NBLK 501 MATDIM = MATDIM + NDIM(IBLK)*NDIM(IBLK) 502 200 CONTINUE 503C 504 AMAX = ZERO 505 N_NAN = 0 506 DO I=1,MATDIM 507 IF ( IS_NAN(AMATRX(I),AMATRX(I)) ) THEN 508 N_NAN = N_NAN + 1 509 ELSE 510 AMAX = MAX( AMAX, ABS(AMATRX(I)) ) 511 END IF 512 END DO 513 IF (N_NAN .GT. 0) WRITE (LUPRI,'(/T6,A,I10,A)') 514 & 'WARNING: matrix contains',N_NAN,' NaN.' 515 IF (AMAX .EQ. ZERO) THEN 516 WRITE (LUPRI,3000) NBLK 517 GO TO 800 518 END IF 519 IF (FFMIN .LE. AMAX .AND. AMAX .LE. FFMAX) THEN 520C use F output format 521 PFMT = '(A1,I7,2X,8F15.8)' 522 ELSE 523C use 1PD output format 524 PFMT = '(A1,I7,2X,1P,8D15.6)' 525 END IF 526C 527 IOFF = 0 528 DO 100 IBLK = 1,NBLK 529 IDIM = NDIM(IBLK) 530 IF (IDIM.EQ.0) GO TO 100 531 I2DIM = IDIM*IDIM 532 DO 10 I=1,I2DIM 533 IF (AMATRX(IOFF+I).NE.ZERO) GO TO 15 534 10 CONTINUE 535 WRITE (LUPRI,1100) IBLK 536 GO TO 90 537 15 CONTINUE 538 WRITE (LUPRI,1200) IBLK 539 LAST = MIN(IDIM,KCOL) 540 DO 2 BEGIN = 1,IDIM,KCOL 541 WRITE (LUPRI,1000) (COLUMN,I,I = BEGIN,LAST) 542 DO 1 K = 1,IDIM 543 IOFFK = IOFF + K 544 IF (MCTL .GE. 10) GO TO 5 ! also print zero rows 545 DO 4 I=BEGIN,LAST 546 IF (AMATRX(IOFFK+(I-1)*IDIM).NE.ZERO) GO TO 5 547 4 CONTINUE 548 GO TO 1 549 5 WRITE (LUPRI,PFMT) CTL,K,(AMATRX(IOFFK+(I-1)*IDIM), 550 * I = BEGIN,LAST) 551 1 CONTINUE 552 2 LAST = MIN(LAST+KCOL,IDIM) 553 90 IOFF = IOFF + I2DIM 554 100 CONTINUE 555C 556 800 CONTINUE 557 WRITE(LUPRI,'(A)') ' ==== End of matrix output ====' 558 RETURN 559 3000 FORMAT (/5X,'All',I3,' blocks zero matrices.') 560 1100 FORMAT (/5X,'*** Block',I3,' zero matrix. ***') 561 1200 FORMAT (/5X,'*** Block',I3,' ***') 562 1000 FORMAT (/10X,8(5X,A6,I4)) 563 END 564C /* Deck output */ 565 SUBROUTINE OUTPUT (AMATRX,ROWLOW,ROWHI,COLLOW,COLHI,ROWDIM,COLDIM, 566 * NCTL,LUPRI) 567C....................................................................... 568C Revised 15-Dec-1983 by Hans Jorgen Aa. Jensen. 569C 16-Jun-1986 hjaaj ( removed Hollerith ) 570C 571C OUTPUT PRINTS A REAL MATRIX IN FORMATTED FORM WITH NUMBERED ROWS 572C AND COLUMNS. THE INPUT IS AS FOLLOWS; 573C 574C AMATRX(',').........MATRIX TO BE OUTPUT 575C 576C ROWLOW..............ROW NUMBER AT WHICH OUTPUT IS TO BEGIN 577C 578C ROWHI...............ROW NUMBER AT WHICH OUTPUT IS TO END 579C 580C COLLOW..............COLUMN NUMBER AT WHICH OUTPUT IS TO BEGIN 581C 582C COLHI...............COLUMN NUMBER AT WHICH OUTPUT IS TO END 583C 584C ROWDIM..............ROW DIMENSION OF AMATRX(',') 585C 586C COLDIM..............COLUMN DIMENSION OF AMATRX(',') 587C 588C NCTL................CARRIAGE CONTROL FLAG; 1 FOR SINGLE SPACE 589C 2 FOR DOUBLE SPACE 590C 3 FOR TRIPLE SPACE 591C hjaaj: negative for 132 col width 592C 593C THE PARAMETERS THAT FOLLOW MATRIX ARE ALL OF TYPE INTEGER. THE 594C PROGRAM IS SET UP TO HANDLE 5 COLUMNS/PAGE WITH A 1P,5D24.15 FORMAT 595C FOR THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, 596C CHANGE FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER 597C OF COLUMNS. 598C 599C AUTHOR; NELSON H.F. BEEBE, QUANTUM THEORY PROJECT, UNIVERSITY OF 600C FLORIDA, GAINESVILLE 601C REVISED; FEBRUARY 26, 1971 602C 603C....................................................................... 604C 605#include "implicit.h" 606 INTEGER ROWLOW,ROWHI,COLLOW,COLHI,ROWDIM,COLDIM,BEGIN,KCOL 607 DIMENSION AMATRX(ROWDIM,COLDIM) 608 CHARACTER*1 ASA(3), BLANK, CTL 609 CHARACTER PFMT*20, COLUMN*8 610 LOGICAL IS_NAN 611 PARAMETER (ZERO=0.D00, KCOLP=5, KCOLN=8) 612 PARAMETER (FFMIN=1.D-3, FFMAX = 1.D3) 613 DATA COLUMN/'Column '/, BLANK/' '/, ASA/' ', '0', '-'/ 614C 615 N_NAN = 0 616 IF (ROWHI.LT.ROWLOW) GO TO 3 617 IF (COLHI.LT.COLLOW) GO TO 3 618C 619 AMAX = ZERO 620 DO 10 J = COLLOW,COLHI 621 DO 10 I = ROWLOW,ROWHI 622 IF ( IS_NAN(AMATRX(I,J),AMATRX(I,J)) ) THEN 623 N_NAN = N_NAN + 1 624 ELSE 625 AMAX = MAX( AMAX, ABS(AMATRX(I,J)) ) 626 END IF 627 10 CONTINUE 628 IF (N_NAN .GT. 0) WRITE (LUPRI,'(/T6,A,I10,A)') 629 & 'WARNING: matrix contains',N_NAN,' NaN.' 630 IF (AMAX .EQ. ZERO) THEN 631 WRITE (LUPRI,'(/T6,A)') 'Zero matrix.' 632 GO TO 3 633 END IF 634 IF (FFMIN .LE. AMAX .AND. AMAX .LE. FFMAX) THEN 635C use F output format 636 PFMT = '(A1,I7,2X,8F15.8)' 637 thrpri = 0.5D-8 638 ELSE 639C use 1PE output format 640 PFMT = '(A1,I7,2X,1P,8E15.6)' 641 thrpri = 1.0D-8*AMAX 642 END IF 643C 644 IF (NCTL .LT. 0) THEN 645 KCOL = KCOLN 646 ELSE 647 KCOL = KCOLP 648 END IF 649 MCTL = ABS(NCTL) 650 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 651 CTL = ASA(MCTL) 652 ELSE 653 CTL = BLANK 654 END IF 655C 656 LAST = MIN(COLHI,COLLOW+KCOL-1) 657 DO 2 BEGIN = COLLOW,COLHI,KCOL 658 WRITE (LUPRI,1000) (COLUMN,I,I = BEGIN,LAST) 659 DO 1 K = ROWLOW,ROWHI 660 IF (MCTL .GE. 10) GO TO 5 ! also print zero rows 661 DO 4 I = BEGIN,LAST 662 IF (abs(AMATRX(K,I)).gt.thrpri) GO TO 5 663 4 CONTINUE 664 GO TO 1 665 5 WRITE (LUPRI,PFMT) CTL,K,(AMATRX(K,I), I = BEGIN,LAST) 666 1 CONTINUE 667 2 LAST = MIN(LAST+KCOL,COLHI) 668 3 WRITE(LUPRI,'(A)') ' ==== End of matrix output ====' 669 RETURN 670 1000 FORMAT (/10X,8(5X,A6,I4)) 671 END 672C /* Deck Ioutput */ 673 SUBROUTINE IOUTPUT(IMATRX,ROWLOW,ROWHI,COLLOW,COLHI,ROWDIM,COLDIM, 674 * NCTL,LUPRI) 675C....................................................................... 676C 04-Jun-2014 by Hans Jorgen Aa. Jensen. 677C Based on OUTPUT. 678C 679C IOUTPUT PRINTS an integer MATRIX IN FORMATTED FORM WITH NUMBERED ROWS 680C AND COLUMNS. THE INPUT IS AS FOLLOWS; 681C 682C IMATRX(',').........integer MATRIX TO BE OUTPUT 683C 684C ROWLOW..............ROW NUMBER AT WHICH OUTPUT IS TO BEGIN 685C 686C ROWHI...............ROW NUMBER AT WHICH OUTPUT IS TO END 687C 688C COLLOW..............COLUMN NUMBER AT WHICH OUTPUT IS TO BEGIN 689C 690C COLHI...............COLUMN NUMBER AT WHICH OUTPUT IS TO END 691C 692C ROWDIM..............ROW DIMENSION OF IMATRX(',') 693C 694C COLDIM..............COLUMN DIMENSION OF IMATRX(',') 695C 696C NCTL................CARRIAGE CONTROL FLAG; 1 FOR SINGLE SPACE 697C 2 FOR DOUBLE SPACE 698C 3 FOR TRIPLE SPACE 699C hjaaj: negative for 132 col width 700C 701C THE PARAMETERS THAT FOLLOW MATRIX ARE ALL OF TYPE INTEGER. THE 702C PROGRAM IS SET UP TO HANDLE 5 COLUMNS/PAGE WITH A 1P,5D24.15 FORMAT 703C FOR THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, 704C CHANGE FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER 705C OF COLUMNS. 706C 707C....................................................................... 708C 709#include "implicit.h" 710 INTEGER ROWLOW,ROWHI,COLLOW,COLHI,ROWDIM,COLDIM,BEGIN,KCOL 711 INTEGER IMATRX(ROWDIM,COLDIM), IM_MAX 712 CHARACTER*1 ASA(3), BLANK, CTL 713 CHARACTER CFMT*20,PFMT*20, COLUMN*8 714 DATA COLUMN/'Column '/, BLANK/' '/, ASA/' ', '0', '-'/ 715C 716 IF (ROWHI.LT.ROWLOW) GO TO 3 717 IF (COLHI.LT.COLLOW) GO TO 3 718C 719 IM_MAX = 0 720 DO 10 J = COLLOW,COLHI 721 DO 10 I = ROWLOW,ROWHI 722 IM_MAX = MAX( IM_MAX, ABS(IMATRX(I,J)) ) 723 10 CONTINUE 724 IF (IM_MAX .EQ. 0) THEN 725 WRITE (LUPRI,'(/T6,A)') 'Zero matrix.' 726 GO TO 3 727 END IF 728 IF (IM_MAX .le. 1000000) THEN 729 CFMT = '(/10X,20(1X,A3,I4))' 730 PFMT = '(A1,I7,2X,20I8)' 731 KCOLP = 8 732 KCOLN = 15 733 ELSE 734 CFMT = '(/10X,8(5X,A6,I4))' 735 PFMT = '(A1,I7,2X,8I15)' 736 KCOLP = 5 737 KCOLN = 8 738 END IF 739C 740 IF (NCTL .LT. 0) THEN 741 KCOL = KCOLN 742 ELSE 743 KCOL = KCOLP 744 END IF 745 MCTL = ABS(NCTL) 746 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 747 CTL = ASA(MCTL) 748 ELSE 749 CTL = BLANK 750 END IF 751C 752 LAST = MIN(COLHI,COLLOW+KCOL-1) 753 DO 2 BEGIN = COLLOW,COLHI,KCOL 754 WRITE (LUPRI,CFMT) (COLUMN,I,I = BEGIN,LAST) 755 DO 1 K = ROWLOW,ROWHI 756 DO 4 I = BEGIN,LAST 757 IF (abs(IMATRX(K,I)).gt.0) GO TO 5 758 4 CONTINUE 759 GO TO 1 760 5 WRITE (LUPRI,PFMT) CTL,K,(IMATRX(K,I), I = BEGIN,LAST) 761 1 CONTINUE 762 2 LAST = MIN(LAST+KCOL,COLHI) 763 3 WRITE(LUPRI,'(A)') ' ==== End of integer matrix output ====' 764 RETURN 765 END 766C /* Deck prtab */ 767 SUBROUTINE PRTAB(NTABLE,TABLE,TEXT,LUPRI) 768C 769C 28-Dec-1989 Hans Joergen Aa. Jensen 770C 771C Purpose: print tables of text (e.g. from input parsing routines) 772C 773 CHARACTER*(*) TABLE(NTABLE), TEXT 774 LTEXT = LEN(TEXT) 775 LTEXT = MIN(70,LTEXT) 776 WRITE (LUPRI,'(//1X,A/1X,70A1/)') TEXT(1:LTEXT),('=',I=1,LTEXT) 777 DO 100 I = 1,NTABLE 778 IF (INDEX(TABLE(I),'XXXX') .EQ. 0) THEN 779 WRITE (LUPRI,'(T6,A)') TABLE(I) 780 END IF 781 100 CONTINUE 782 WRITE (LUPRI,'()') 783 RETURN 784 END 785C /* Deck prvec */ 786 SUBROUTINE PRVEC(NDIM,VEC,INCVEC,THRESH,MAXLIN,LUOUT) 787C 788C 19-Aug-1989 Hans Joergen Aa. Jensen 789C 790C print VEC(1:NDIM*INCVEC:INCVEC) 791C 792C NDIM : Number of elements in vector VEC 793C VEC(:) : Vector to be printed 794C INCVEC : Increment between each element in vector VEC 795C THRESH : Print threshold for vector with unit norm 796C (if THRESH .lt. 0 then -THRESH is used without 797C renormalization). 798C MAXLIN : max. lines of output with vector elements 799C LUOUT : output unit 800C 801C 802#include "implicit.h" 803 DIMENSION VEC(*) 804 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 805 PARAMETER (D1LOW = D1 - 1.0D-10, D1HGH = D1 + 1.0D-10) 806 LOGICAL VSCALE 807C 808C Test input 809C 810 NERR = 0 811 IF (NDIM .LE. 0) THEN 812 WRITE (LUOUT,'(/5X,A)') 813 & 'No print from PRVEC because NDIM .le. 0' 814 NERR = NERR + 1 815 END IF 816 IF (INCVEC .LE. 0) THEN 817 WRITE (LUOUT,'(/5X,A)') 818 & 'No print from PRVEC because INCVEC .le. 0' 819 NERR = NERR + 1 820 END IF 821 IF (NERR .GT. 0) RETURN 822C 823C 824C 825 C2NRM = DNRM2(NDIM,VEC,INCVEC) 826 IF (THRESH .LE. D0) THEN 827 THRES2 = -THRESH 828 VSCALE = .FALSE. 829 ELSE 830 THRES2 = THRESH * C2NRM 831 VSCALE = (C2NRM .LT. D1LOW .OR. C2NRM .GT. D1HGH) 832 END IF 833 IF (VSCALE) THEN 834 WRITE (LUOUT,'(/2A,1P,D12.4)') 835 & ' Print of vector elements (vector scaled to unit norm)', 836 & ' larger than',ABS(THRESH) 837 SCALE = D1 / C2NRM 838 ELSE 839 WRITE (LUOUT,'(/A,1P,D10.2)') 840 & ' Print of vector elements larger than',ABS(THRESH) 841 SCALE = D1 842 END IF 843C 844 IPR = 0 845 IZER = 0 846 C2SUM = D0 847C 848 DO 300 I = 1, NDIM 849 NA = 1 + (I-1)*INCVEC 850 IF (ABS(VEC(NA)).LE.THRES2 .OR. IPR .GE. MAXLIN) THEN 851 C2SUM = C2SUM + VEC(NA)*VEC(NA) 852 IF (VEC(NA) .EQ. D0) IZER = IZER + 1 853 ELSE 854 IF (MOD(IPR,5) .EQ. 0) WRITE (LUOUT,'()') 855 IPR = IPR + 1 856 IF (VSCALE) THEN 857 WRITE(LUOUT,50)I,VEC(NA),SCALE*VEC(NA) 858 ELSE 859 WRITE(LUOUT,60)I,VEC(NA) 860 END IF 861 50 FORMAT(3X,'Element',I10,3X,'coefficient',1P,D10.2,0P, 862 & 3X,'scaled to unit norm',F10.6) 863 60 FORMAT(6X,'Element',I12,3X,'coefficient',F20.10) 864 END IF 865 300 CONTINUE 866 IF (IPR .GE. MAXLIN) THEN 867C 868C *** We have reached the print limit 869C 870 WRITE (LUOUT,910) IPR 871 END IF 872 910 FORMAT(/'INFO: Print limit of',I6,' elements has been reached.') 873 C2SUM = SQRT(C2SUM) 874 IF (IZER .EQ. NDIM) THEN 875 WRITE (LUOUT,920) NDIM 876 ELSE 877 WRITE (LUOUT,930) NDIM,NDIM-IPR,IZER,C2SUM,C2NRM 878 END IF 879 920 FORMAT(/' Length of vector :',I10, 880 * /' All elements are zero.',/) 881 930 FORMAT(/' Length of vector :',I10, 882 * /' Number of elements not printed :',I10, 883 * /' Number of zero elements :',I10, 884 * /' Total norm of coefficients not printed:',F10.6, 885 * /' (the coefficients are normalized to ',F10.6,')',/) 886 RETURN 887C 888C End of PRVEC 889C 890 END 891C /* Deck prmgn */ 892 SUBROUTINE PRMGN(NDIM,VEC,INCVEC,NPOT,LUOUT) 893C 894C 19-Aug-1989 hjaaj 895C 896C NDIM : Number of elements in vector VEC 897C VEC(:) : Vector to be printed 898C INCVEC : Increment between each element in vector VEC 899C NPOT : IBASE**-NPOT is lowest magnitude considered 900C LUOUT : output unit 901C 902C 903#include "implicit.h" 904 DIMENSION VEC(*) 905 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 906 PARAMETER (IBASE = 10) 907C IBASE : Base for magnitude 908C 909 FACMGN = IBASE 910 FACMGN = D1 / FACMGN 911C 912C Test input 913C 914 NERR = 0 915 IF (NDIM .LE. 0) THEN 916 WRITE (LUOUT,*) 'No print from PRMGN because NDIM .le. 0' 917 NERR = NERR + 1 918 END IF 919 IF (INCVEC .LE. 0) THEN 920 WRITE (LUOUT,*) 'No print from PRMGN because INCVEC .le. 0' 921 NERR = NERR + 1 922 END IF 923 IF (NERR .GT. 0) RETURN 924C 925C 926C 927 C2NRM = DNRM2(NDIM,VEC,INCVEC) 928 IZER = 0 929 NLAST = NDIM*INCVEC 930 VECMAX = D0 931 DO 100 NA = 1, NLAST, INCVEC 932 IF (VEC(NA) .EQ. D0) THEN 933 IZER = IZER + 1 934 ELSE 935 VECMAX = MAX(VECMAX, ABS(VEC(NA)) ) 936 END IF 937 100 CONTINUE 938C 939C.... SIZE OF COEFFICIENTS 940C 941 WRITE(LUOUT,'(/A/A//A,2I12)') 942 & ' Magnitude of coefficients ', 943 & ' ===========================', 944 & ' Number of elements and increment:',NDIM,INCVEC 945 IF (IZER .EQ. NDIM) THEN 946 WRITE(LUOUT,'(/A)') ' All elements are zero.' 947 GO TO 9999 948 END IF 949 WRITE(LUOUT,'(/A,1P,2E23.12)') 950 & ' Maximum abs value and vector norm:',VECMAX,C2NRM 951C 952 IF (C2NRM .GT. (D1 + 1.0 D-10) ) THEN 953 WRITE (LUOUT,'(/A,1P,E23.12,A,/A)') 954 & ' (Vector norm is',C2NRM,')', 955 & ' (Range is scaled to a vector norm of 1)' 956 END IF 957 XMIN = MAX(C2NRM,D1) 958C ... max possible coefficient is C2NRM 959C 960 WRITE (LUOUT,'(/2A,/2A)') 961 & ' Range # of elements', 962 & ' Norm squared Accumulated', 963 & ' ------------- --------------', 964 & ' --------------- ---------------' 965C 966C Output will look like: 967C 968C Range # of elements Norm squared Accumulated 969C ------------- -------------- --------------- --------------- 970C 10- 1 to 1 2 1.12345678E-01 1.12345678E-01 971C 972 C2NRM = D0 973 IDET = IZER 974C. LOOP OVER INTERVALS 975 IPOT = 0 976 200 CONTINUE 977 CLNORM = D0 978 INRANG = 0 979 XMAX = XMIN 980 IF (IPOT .EQ. 0) XMAX = 2*XMAX 981C ... to make certain we don't omit any because of round-off 982 XMIN = XMIN * FACMGN 983 DO 300 NA = 1, NLAST, INCVEC 984 IF( ABS(VEC(NA)) .LE. XMAX .AND. 985 & ABS(VEC(NA)) .GT. XMIN ) THEN 986 INRANG = INRANG + 1 987 CLNORM = CLNORM + VEC(NA) ** 2 988 END IF 989 300 CONTINUE 990 C2NRM = C2NRM + CLNORM 991C 992C 993 IF (INRANG .GT. 0) THEN 994 IF (IPOT .EQ. 0) THEN 995 WRITE(LUOUT,'(I3,A,I14,1P,2E20.8)') 996 & IBASE,'- 1 to 1 ',INRANG,CLNORM,C2NRM 997 ELSE 998 WRITE(LUOUT,'(I3,A,I2,A,I3,A,I2,I14,1P,2E20.8)') 999 & IBASE,'-',IPOT+1,' to',IBASE,'-',IPOT, 1000 & INRANG,CLNORM,C2NRM 1001 END IF 1002 END IF 1003C 1004C 1005 IDET = IDET + INRANG 1006 IPOT = IPOT + 1 1007 IF( IDET .LT. NDIM .AND. IPOT .LT. NPOT ) GOTO 200 1008C 1009 ISML = NDIM - IDET 1010 IF (ISML .GT. 0) WRITE (LUOUT,'(A,I13)') ' Other non-zero ',ISML 1011 IF (IZER .GT. 0) WRITE (LUOUT,'(A,I13)') ' Exact zero ',IZER 1012 WRITE (LUOUT,'(/A,I13/)') ' Total ',NDIM 1013C 1014C 1015 9999 CONTINUE 1016 RETURN 1017C 1018C End of PRMGN 1019C 1020 END 1021C /* Deck pnzvec */ 1022 FUNCTION PNZVEC(NDIM,VEC,INCVEC,THRZER) 1023C 1024C Copyright 24-Mar-1993 Hans Joergen Aagaard Jensen 1025C return % non-zero elements of VEC(1:NDIM*INCVEC:INCVEC) 1026C 1027C NDIM : Number of elements in vector VEC 1028C VEC(:) : Vector 1029C INCVEC : Increment between each element in vector VEC 1030C THRZER : Threshold for zero elements 1031C 1032C 1033#include "implicit.h" 1034 DIMENSION VEC(*) 1035#include "priunit.h" 1036 PARAMETER (D0 = 0.0D0, D100 = 100.0D0) 1037C 1038C Test input 1039C 1040 NERR = 0 1041 IF (NDIM .LE. 0) THEN 1042 WRITE (LUPRI,'(/5X,A,I10)') 1043 & 'No percentage from PNZVEC because NDIM .le. 0; NDIM =',NDIM 1044 NERR = NERR + 1 1045 END IF 1046 IF (INCVEC .LE. 0) THEN 1047 WRITE (LUPRI,'(/5X,A,I10)') 1048 & 'No percentage from PNZVEC because INCVEC .le. 0; INCVEC =', 1049 & INCVEC 1050 NERR = NERR + 1 1051 END IF 1052 IF (NERR .GT. 0) THEN 1053 PNZVEC = -D100 1054 RETURN 1055 END IF 1056C 1057 THRESH = MAX(D0,THRZER) 1058 NNZ = 0 1059 DO 100 I = 1,NDIM*INCVEC,INCVEC 1060 IF (ABS(VEC(I)) .GT. THRESH) NNZ = NNZ + 1 1061 100 CONTINUE 1062 DNZ = NNZ 1063 DNDIM = NDIM 1064 PNZVEC = DNZ*D100 / DNDIM 1065 RETURN 1066 END 1067C /* Deck coutput */ 1068 SUBROUTINE COUTPUT(AMATRX,ROWLOW,ROWHI,COLLOW,COLHI,ROWDIM,COLDIM, 1069 * NCTL,LUPRI) 1070C....................................................................... 1071C Revised 15-Dec-1983 by Hans Jorgen Aa. Jensen. 1072C 16-Jun-1986 hjaaj ( removed Hollerith ) 1073C 20-Mar-2001 ov complex matrix verison 1074C 1075C OUTPUT PRINTS A REAL MATRIX IN FORMATTED FORM WITH NUMBERED ROWS 1076C AND COLUMNS. THE INPUT IS AS FOLLOWS; 1077C 1078C AMATRX(',').........MATRIX TO BE OUTPUT 1079C 1080C ROWLOW..............ROW NUMBER AT WHICH OUTPUT IS TO BEGIN 1081C 1082C ROWHI...............ROW NUMBER AT WHICH OUTPUT IS TO END 1083C 1084C COLLOW..............COLUMN NUMBER AT WHICH OUTPUT IS TO BEGIN 1085C 1086C COLHI...............COLUMN NUMBER AT WHICH OUTPUT IS TO END 1087C 1088C ROWDIM..............ROW DIMENSION OF AMATRX(',') 1089C 1090C COLDIM..............COLUMN DIMENSION OF AMATRX(',') 1091C 1092C NCTL................CARRIAGE CONTROL FLAG; 1 FOR SINGLE SPACE 1093C 2 FOR DOUBLE SPACE 1094C 3 FOR TRIPLE SPACE 1095C hjaaj: negative for 132 col width 1096C 1097C THE PARAMETERS THAT FOLLOW MATRIX ARE ALL OF TYPE INTEGER. THE 1098C PROGRAM IS SET UP TO HANDLE 5 COLUMNS/PAGE WITH A 1P,5D24.15 FORMAT 1099C FOR THE COLUMNS. IF A DIFFERENT NUMBER OF COLUMNS IS REQUIRED, 1100C CHANGE FORMATS 1000 AND 2000, AND INITIALIZE KCOL WITH THE NEW NUMBER 1101C OF COLUMNS. 1102C 1103C AUTHOR; NELSON H.F. BEEBE, QUANTUM THEORY PROJECT, UNIVERSITY OF 1104C FLORIDA, GAINESVILLE 1105C REVISED; FEBRUARY 26, 1971 1106C 1107C....................................................................... 1108C 1109 IMPLICIT DOUBLE COMPLEX (A-H,O-Z) 1110 INTEGER ROWLOW,ROWHI,COLLOW,COLHI,ROWDIM,COLDIM,BEGIN,KCOL 1111 DIMENSION AMATRX(ROWDIM,COLDIM) 1112 COMPLEX C4_test ! hjaaj: for testing elements for print 1113 REAL*4 R4_test,R4_thrpri ! because no double complex absolute value :-( 1114 CHARACTER*1 ASA(3), BLANK, CTL 1115 CHARACTER PFMT*30, COLUMN*8 1116 DOUBLE PRECISION ZERO, FFMIN, FFMAX 1117 PARAMETER (ZERO=0.D00, KCOLP=4, KCOLN=6) 1118 PARAMETER (FFMIN=1.D-3, FFMAX = 1.D3) 1119 DATA COLUMN/'Column '/, BLANK/' '/, ASA/' ', '0', '-'/ 1120 DOUBLE PRECISION AMAX 1121C 1122 IF (ROWHI.LT.ROWLOW) GO TO 3 1123 IF (COLHI.LT.COLLOW) GO TO 3 1124C 1125 AMAX = ZERO 1126 DO 10 J = COLLOW,COLHI 1127 DO 10 I = ROWLOW,ROWHI 1128 AMAX = MAX( AMAX, ABS(AMATRX(I,J)) ) 1129 10 CONTINUE 1130 IF (AMAX .EQ. ZERO) THEN 1131 WRITE (LUPRI,'(/T6,A)') 'Zero matrix.' 1132 GO TO 3 1133 END IF 1134 IF (FFMIN .LE. AMAX .AND. AMAX .LE. FFMAX) THEN 1135C use F output format 1136 PFMT = '(A1,I7,2X,6(F17.8,F13.8))' 1137 thrpri = 0.5D-6 1138 ELSE 1139C use 1PD output format 1140 PFMT = '(A1,I7,2X,1P,6(D17.6,D13.6))' 1141 thrpri = 1.0D-6*AMAX 1142 END IF 1143C 1144 IF (NCTL .LT. 0) THEN 1145 KCOL = KCOLN 1146 ELSE 1147 KCOL = KCOLP 1148 END IF 1149 MCTL = ABS(NCTL) 1150 IF ((MCTL.LE.3).AND.(MCTL.GT.0)) THEN 1151 CTL = ASA(MCTL) 1152 ELSE 1153 CTL = BLANK 1154 END IF 1155C 1156 R4_thrpri = thrpri 1157 LAST = MIN(COLHI,COLLOW+KCOL-1) 1158 DO 2 BEGIN = COLLOW,COLHI,KCOL 1159 WRITE (LUPRI,1000) (COLUMN,I,I = BEGIN,LAST) 1160 DO 1 K = ROWLOW,ROWHI 1161 IF (MCTL .GE. 10) GO TO 5 ! also print zero rows 1162 DO 4 I = BEGIN,LAST 1163 C4_test = AMATRX(K,I) ! hjaaj: all this because no intrinsic routine to 1164 R4_test = cabs(C4_test) ! calculate absolute value of double complex number 1165 IF (R4_test .gt. R4_thrpri) GO TO 5 1166 4 CONTINUE 1167 GO TO 1 1168 5 WRITE (LUPRI,PFMT) CTL,K,(AMATRX(K,I), I = BEGIN,LAST) 1169 1 CONTINUE 1170 2 LAST = MIN(LAST+KCOL,COLHI) 1171 3 WRITE(LUPRI,'(A)') ' ==== End of matrix output ====' 1172 RETURN 1173 1000 FORMAT (/10X,6(' real < ',A6,I4,' > imag')) 1174 END 1175C === Deck is_nan =============== 1176 LOGICAL FUNCTION IS_NAN(XA,XB) 1177C 1178C May 2010, Hans Joergen Aa. Jensen 1179C Purpose: IS_NAN(X,X) is true iff X is NAN 1180C 1181 REAL*8 XA, XB 1182 IS_NAN = XA .NE. XB 1183 RETURN 1184 END 1185C === Deck n_nan =============== 1186 INTEGER FUNCTION N_NAN(N,AVEC) 1187C 1188C May 2010, Hans Joergen Aa. Jensen 1189C Purpose: N_NAN returns number of NAN entries in AVEC(1:N) 1190C 1191 LOGICAL IS_NAN 1192 INTEGER N, I, J 1193 REAL*8 AVEC(N) 1194 J = 0 1195 DO I = 1,N 1196 IF ( IS_NAN(AVEC(I),AVEC(I)) ) J = J + 1 1197 END DO 1198 N_NAN = J 1199 RETURN 1200 END 1201C === Deck n_not_equal =============== 1202 INTEGER FUNCTION N_NOT_EQUAL(N,AVEC,BVEC) 1203C 1204C May 2010, Hans Joergen Aa. Jensen 1205C Purpose: Number of not equal elements in AVEC and BVEC 1206C 1207C NOTE: N_NOT_EQUAL(N,DA,DA) is number of NAN elements in DA(1:N) 1208C 1209 INTEGER N, I, J 1210 REAL*8 AVEC(N), BVEC(N) 1211 J = 0 1212 DO I = 1,N 1213 IF ( AVEC(I) .NE. BVEC(I) ) J = J + 1 1214 END DO 1215 N_NOT_EQUAL = J 1216 RETURN 1217 END 1218C -------- end of printpkg.F ----------- 1219