1* 2* $Id$ 3* 4c 5c Set of routines to cache AO integrals for the 4-index 6c module using exclusive-access local disk storage. 7c Integrals are retrieved in subsequent passes through 8c AO list. 9c 10c Maintain state via two common block variables. 11c 12c moao_ipass < 0 (default) disk caching disabled 13c = 0 saving to disk, initial pass 14c > 0 retrieving from disk 15c 16c moao_fd < 0 io not initialized 17c >= 0 io is happening 18c 19c Unless explicitly initialized, there will be no disk 20c caching. 21c 22c Since exclusive-access files are used, the same subset of 23c AO integrals that was initially generated on a processor 24c will be retrieved in subsequent passes. Thus, the task scheduling 25c from the initial run must be emulated. 26c Use a wrapper around nxtask() and record the task number on 27c each disk record. 28c 29c 30c 31c 32c 33c Following routines use the common block 34c 35 block data moints_moao_block 36 implicit none 37#include "cmointsmoao.fh" 38 data moao_ipass/-1/ 39 data moao_fd/-1/ 40 data moao_lbuf/-1/ 41 data moao_fname/' '/ 42 end 43 44 45 46 47 48c 49c Wrapper around nxtask so we can save task numbers 50c with the IO records and emulate task scheduling 51c when retrieved 52c 53 integer function moints_nxttask( numnodes ) 54 implicit none 55#include "tcgmsg.fh" 56#include "global.fh" 57#include "cmointsmoao.fh" 58 integer numnodes 59 integer nxtask 60 external nxtask 61 62 if (moao_ipass.gt.0) then 63 moints_nxttask = moao_tasknum 64 if (numnodes.lt.0) call ga_sync() 65 else 66 moints_nxttask = nxtask( numnodes, 1 ) 67 if (moao_ipass.eq.0) then 68 moao_tasknum = moints_nxttask 69 endif 70 endif 71 return 72 end 73 74 75 76c 77c Return complete SSBB block of integrals if cached on disk 78c 79 logical function moints_gblk_fromdisk( blkid, ish, jsh, 80 $ kshlo, lshlo, 81 $ ilo, ihi, jlo, jhi, 82 $ kblo, kbhi, lblo, lbhi, 83 $ ssbb ) 84 implicit none 85#include "errquit.fh" 86 integer blkid 87 integer ish, jsh 88 integer kshlo 89 integer lshlo 90 integer ilo, ihi 91 integer jlo, jhi 92 integer kblo, kbhi 93 integer lblo, lbhi 94 double precision ssbb( lblo:lbhi, kblo:kbhi, jlo:jhi, ilo:ihi ) 95c 96#include "mafdecls.fh" 97#include "cmointsmoao.fh" 98c 99#ifdef OLD_AODISK 100 character*8 buidstr 101 integer recnum 102 logical moints_iorec_next 103 external moints_iorec_next 104#endif 105 integer ssbblen 106c$$$ DOUBLE PRECISION DABSSUM 107c$$$ EXTERNAL DABSSUM 108 109 moints_gblk_fromdisk = .false. 110 if (moao_ipass.lt.0) return 111 if ((moao_fd.ge.0).and.(moao_ipass.gt.0)) then 112 113#ifdef OLD_AODISK 114 recnum = 0 115 ssbblen = (ihi-ilo+1)*(jhi-jlo+1)*(kbhi-kblo+1)*(lbhi-lblo+1) 116 call dfill( ssbblen, 0.d0, ssbb, 1 ) 117 do while (moints_iorec_next( blkid, buidstr )) 118 if (moao_issparse.eq.1) then 119 call moints_aodisk_iorec2sprs_old( ilo, ihi, jlo, jhi, 120 $ kblo, kbhi, lblo, lbhi, 121 $ ssbb, moao_reclen, moao_lwidth, 122 $ dbl_mb(moao_klabrec), 123 $ dbl_mb(moao_kvalrec) ) 124 moints_gblk_fromdisk = .true. 125 else 126 call errquit('moints_gblk_disk: dense write not ready',0, 127 & INT_ERR) 128 endif 129 recnum = recnum + 1 130 enddo 131#else 132 ssbblen = (ihi-ilo+1)*(jhi-jlo+1)*(kbhi-kblo+1)*(lbhi-lblo+1) 133 call dfill( ssbblen, 0.d0, ssbb, 1 ) 134 call moints_aodisk_iorec2sprs( moao_fd, moao_fptr, 135 $ moao_tasknum, 136 $ ilo, ihi, jlo, jhi, 137 $ kblo, kbhi, lblo, lbhi, 138 $ ssbb, moao_hdrp, moao_buflen, 139 $ dbl_mb(moao_kbuf), moao_nrec ) 140 moints_gblk_fromdisk = .true. 141#endif 142C MOAO_CUMUL = MOAO_CUMUL + DABSSUM( SSBBLEN, SSBB ) 143c$$$ WRITE(6,771) ISH, JSH, KSHLO, LSHLO, MOAO_TASKNUM, GA_NODEID(), 144c$$$ $ BLKID, DABSSUM(SSBBLEN,SSBB) 145c$$$ 771 FORMAT('qqqq-',4I4,3X,I5,I3,5X,I6,5X,F20.6) 146 endif 147 return 148 end 149 150 151 152 153c 154c Save SSBB block to disk if caching enabled. 155c 156 logical function moints_gblk_todisk( blkid, ish, jsh, 157 $ kshlo, lshlo, 158 $ ilo, ihi, jlo, jhi, 159 $ kblo, kbhi, lblo, lbhi, 160 $ ssbb ) 161 implicit none 162#include "errquit.fh" 163 integer blkid 164 integer ish, jsh 165 integer kshlo, lshlo 166 integer ilo, ihi 167 integer jlo, jhi 168 integer kblo, kbhi 169 integer lblo, lbhi 170 double precision ssbb( lblo:lbhi, kblo:kbhi, jlo:jhi, ilo:ihi ) 171c 172#include "mafdecls.fh" 173#include "cmointsmoao.fh" 174c 175c 176 logical sparse ! for the moment only sparse case 177 data sparse/.true./ 178c 179c$$$ INTEGER SSBBLEN 180c$$$ DOUBLE PRECISION DABSSUM 181c$$$ EXTERNAL DABSSUM 182 183 moints_gblk_todisk = .false. 184 if (moao_ipass.lt.0) return 185 if ((moao_fd.ge.0).and.(moao_ipass.eq.0)) then 186 if (sparse) then 187#ifdef OLD_AODISK 188 call moints_aodisk_sprs2iorec_old( moao_fd, moao_fptr, 189 $ moao_tasknum, blkid, 190 $ ilo, ihi, jlo, jhi, 191 $ kblo, kbhi, lblo, lbhi, 192 $ ssbb, moao_spreclen, moao_lwidth, 193 $ dbl_mb(moao_klabrec), 194 $ dbl_mb(moao_kvalrec), 195 $ moao_buflen, dbl_mb(moao_kbuf) ) 196#else 197 call moints_aodisk_sprs2iorec( moao_fd, moao_fptr, 198 $ moao_tasknum, blkid, 199 $ ilo, ihi, jlo, jhi, 200 $ kblo, kbhi, lblo, lbhi, 201 $ ssbb, moao_recptr, moao_buflen, 202 $ dbl_mb(moao_kbuf), moao_nrec ) 203#endif 204 else 205 call errquit('moints_gblk_todisk: dense read not ready',0, 206 & INT_ERR) 207 endif 208 moints_gblk_todisk = .true. 209c 210C MOAO_CUMUL = MOAO_CUMUL + DABSSUM( SSBBLEN, SSBB ) 211c$$$ SSBBLEN = (IHI-ILO+1)*(JHI-JLO+1)*(KBHI-KBLO+1)*(LBHI-LBLO+1) 212c$$$ WRITE(6,771) ISH, JSH, KSHLO, LSHLO, MOAO_TASKNUM, GA_NODEID(), 213c$$$ $ BLKID, DABSSUM(SSBBLEN,SSBB) 214c$$$ 771 FORMAT('qqqq+',4I4,3X,I5,I3,5X,I6,5X,F20.6) 215c 216 endif 217 return 218 end 219 220 221 222 223 224 225c 226c Called by application 227c Enable disk caching & check for existing file and open file 228c 229 logical function moints_aodisk_init( odisk, oreuse ) 230 implicit none 231#include "errquit.fh" 232#include "mafdecls.fh" 233#include "global.fh" 234#include "eaf.fh" 235#include "mointsmoaodef.fh" 236#include "cmointsmoao.fh" 237 logical odisk ! [input] toggle caching 238 logical oreuse ! [input] toggle reuse of existing file 239c 240 logical fexist 241 integer stat 242c 243c first time through --- enable caching TO disk 244c otherwise must be caching FROM disk 245c 246 moints_aodisk_init = .false. 247 if (.not.(odisk)) then 248 moao_ipass = -1 249 moints_aodisk_init = .true. 250 return 251 endif 252 if (moao_ipass.lt.0) moao_ipass = 0 253c 254c if AO file exists AND reuse enabled then set ipass > 0 255c 256 if (moao_fname.eq.' ') then 257 call util_file_name( 'moao', .true., .true., moao_fname ) 258 endif 259c 260 inquire(file=moao_fname, exist=fexist) ! Equivalent is eaf_stat ? 261 if (oreuse) then 262 stat=0 263 if (fexist) stat=1 264! igop('*') detects the presence of one or more stat=0 inputs 265#if NWCHEM_USE_IGOP_PROD 266 call ga_igop( 481, stat, 1, '*') 267#else 268 call ga_igop( 481, stat, 1, 'min') 269#endif 270 if (stat.eq.1) then 271 moao_ipass = 1 272 if (ga_nodeid().eq.0) write(6,331) 273 331 format(10x,'Existing AO integral file will be reused') 274 else if (fexist) then 275 call util_file_unlink(moao_fname) ! Only some nodes had it 276 endif 277 else if (fexist) then 278 call util_file_unlink(moao_fname) ! No reuse but file there. 279 endif 280c 281c open AO integral disk cache 282c 283 if (moao_fd.lt.0) then 284 stat = eaf_open( moao_fname, EAF_RW, moao_fd ) 285 if (stat.ne.0) call errquit( 286 $ 'moints_aodisk_init: cannot open ao file',0, DISK_ERR) 287 endif 288 moints_aodisk_init = .true. 289 return 290 end 291 292 293 294 295 296 297 298c 299c This is called internally by moints 300c Initialize disk caching for both saving and retrieving 301c Allocate buffer space and read in first record 302c 303 logical function moints_aodisk_prep( ) 304 implicit none 305#include "errquit.fh" 306#include "mafdecls.fh" 307#include "eaf.fh" 308#include "mointsmoaodef.fh" 309#include "cmointsmoao.fh" 310 integer stat 311 integer blkinfo(MOINTS_NBLKINFO) 312 integer vlsize, dhdrsize, vallen, bufbytes 313 integer nblk,nexpand 314 double precision sparsesize, densesize 315 common/aodisk_stat/nblk,nexpand,sparsesize,densesize 316 317 moints_aodisk_prep = .false. 318 if (moao_ipass.lt.0) then 319 moints_aodisk_prep = .true. 320 return 321 endif 322 moao_nrec = 1 323c 324c allocate io buffer 325c header comes first 326c values and labels start after offset of header length 327c 328 if (moao_lbuf.lt.0) then 329 moao_lwidth = 8/ma_sizeof(MT_INT, 1, MT_BYTE) 330 vlsize = ma_sizeof(MT_INT, moao_lwidth, MT_BYTE) + 331 $ ma_sizeof(MT_DBL, 1, MT_BYTE) ! width of value + labels 332 dhdrsize = ma_sizeof(MT_INT, MOINTS_NBLKINFO, MT_DBL ) 333 if (.not.ma_alloc_get(MT_DBL, MOINTS_IOBUFLEN, 'moints iobuf', 334 $ moao_lbuf, moao_kbuf)) call errquit( 335 $ 'moints_aodisk_prep: cannot allocate io buffer',0, 336 & DISK_ERR) 337 moao_buflen = MOINTS_IOBUFLEN ! bufflen (double words) 338 vallen = moao_buflen - dhdrsize ! value len -= hdr length 339 moao_spreclen = vallen/(vlsize/ma_sizeof(MT_DBL,1,MT_BYTE)) ! leng 340 moao_kvalrec = moao_kbuf + dhdrsize ! values after header 341 moao_klabrec = moao_kbuf + dhdrsize + moao_spreclen ! labels after values 342 endif 343#ifndef OLD_AODISK 344c 345c Initialize first IO record 346c 347 if (moao_ipass.eq.0) then 348 moao_recptr = 2 349 call icopy(1,moao_recptr,1,dbl_mb(moao_kbuf),1) 350 endif 351#endif 352c 353c rewind AO integral cache unit 354c 355 moao_fptr = 1.d0 356 moao_eof = .false. 357c 358c if READ mode -- retrieve first record 359c 360 if (moao_ipass.gt.0) then 361 bufbytes = ma_sizeof(MT_DBL,moao_buflen,MT_BYTE) 362 stat = eaf_read( moao_fd, moao_fptr, dbl_mb(moao_kbuf), 363 $ bufbytes ) 364 moao_fptr = moao_fptr + bufbytes 365#ifdef OLD_AODISK 366 if (stat.eq.0) then 367 call icopy(MOINTS_NBLKINFO, dbl_mb(moao_kbuf), 1, blkinfo, 1) 368 moao_touch = .false. 369 moao_tasknum = blkinfo(1) 370 moao_reclen = blkinfo(2) 371 moao_issparse = blkinfo(3) 372 moao_blkid = blkinfo(4) 373c$$$ WRITE(6,324) moao_tasknum, moao_reclen, moao_blkid 374c$$$ 324 FORMAT(' First record: tasknum=',I5,' RecLen=',I5, 375c$$$ $ ' Blkid=',I5) 376 else 377 moao_ipass = 0 ! this is peculiar but not an error, just reset 378 endif 379#else 380 if (stat.eq.0) then 381 call icopy(MOINTS_NBLKINFO,dbl_mb(moao_kbuf+1),1,blkinfo,1) 382 moao_touch = .false. 383 moao_tasknum = blkinfo(4) 384 moao_blkid = blkinfo(5) 385 moao_issparse = blkinfo(6) 386 moao_recptr = 2 387 moao_hdrp = 2 388 else 389 moao_ipass = 0 390 endif 391#endif 392 endif 393 moints_aodisk_prep = .true. 394 MOAO_CUMUL = 0.d0 395 NEXPAND = 0 396 NBLK = 0 397 SPARSESIZE = 0.d0 398 DENSESIZE = 0.d0 399 return 400 end 401 402 403 404 405c 406c Called internally by moints between passes 407c Free IO buffer and increment pass count 408c Keep NO state information except moao_ipass and moao_fd 409c between passes. 410c 411 subroutine moints_aodisk_tidy() 412 implicit none 413#include "errquit.fh" 414#include "mafdecls.fh" 415#include "global.fh" 416#include "util.fh" 417#include "mointsmoaodef.fh" 418#include "cmointsmoao.fh" 419 double precision fcompress 420 integer nblk, nexpand 421 double precision sparsesize, densesize 422 common/aodisk_stat/nblk,nexpand,sparsesize,densesize 423 424 if (moao_ipass.lt.0) return 425#ifndef OLD_AODISK 426 if (moao_ipass.eq.0) then 427 call moints_aodisk_flushiorec( moao_fd, moao_fptr, ! flush last record 428 $ moao_buflen, dbl_mb(moao_kbuf) ) 429 430 fcompress = (sparsesize/densesize)*100.d0 431 if (util_print('compress stats',print_high).and. 432 $ ga_nodeid().eq.0) then 433 write(6,911) nblk, nexpand, fcompress 434 911 format(' Total blocks = ',I8,/, 435 $ ' Expanded blocks = ',I8,/, 436 $ ' Compression factor = ',F8.1,'%') 437 endif 438 439c$$$ CALL MOINTS_AODISK_CHECKFILE( MOAO_FD, MOAO_BUFLEN, ! verify file for debugging 440c$$$ $ DBL_MB(MOAO_KBUF) ) 441 endif 442#endif 443 if (moao_lbuf.ge.0) then 444 if (.not.ma_free_heap(moao_lbuf)) call errquit( 445 $ 'moints_closeaodisk: cannot free io buffer',0, MA_ERR) 446 moao_lbuf = -1 447 endif 448 moao_ipass = moao_ipass + 1 449 450c$$$ WRITE(6,991) MOAO_IPASS, MOAO_CUMUL, MOAO_NREC 451c$$$ 991 format(' File pass=',i3,/, 452c$$$ $ ' Cumulative value:',f20.6, 453c$$$ $ ' Records read:',i5) 454 455 return 456 end 457 458 459 460 461c 462c Called by application 463c Close I/O unit with option to save 464c Reset ipass count 465c 466 subroutine moints_aodisk_close( osave ) 467 implicit none 468#include "errquit.fh" 469#include "mafdecls.fh" 470#include "eaf.fh" 471#include "util.fh" 472#include "mointsmoaodef.fh" 473#include "cmointsmoao.fh" 474 logical osave 475c 476 integer stat 477 478 if (moao_ipass.lt.0) return 479 moao_ipass = 0 480 if (moao_fd.ge.0) then 481 if (util_print('ao disk stats',print_high)) 482 $ call eaf_print_stats(moao_fd) 483 stat = eaf_close(moao_fd) 484 if (stat.ne.0) call errquit( 485 $ 'moints_closeaodisk: cannot close file',0, DISK_ERR) 486 if (.not.(osave)) then 487 stat = eaf_delete(moao_fname) 488 if (stat.ne.0) call errquit( 489 $ 'moints_closeaodisk: cannot delete file',0, DISK_ERR) 490 endif 491 moao_fd = -1 492 endif 493 494 return 495 end 496 497 498 499 500 501 502 503 504c 505c ====================================================================== 506c 507c 508c Utility routines 509c 510c 511c ====================================================================== 512c 513c 514#ifdef OLD_AODISK 515c 516c 517c Pack the dense SSBB block into sparse form 518c with 16 bits per label 519c 520c 521 subroutine moints_aodisk_sprs2iorec_old( fd, fptr, tasknum, blkid, 522 $ ilo, ihi, jlo, jhi, 523 $ kblo, kbhi, lblo, lbhi, 524 $ ssbb, reclen, lwidth, 525 $ iolab, ioval, 526 $ iobuflen, iobuf ) 527 implicit none 528 integer fd 529 double precision fptr 530 integer tasknum 531 integer blkid 532 integer ilo, ihi, jlo, jhi 533 integer kblo, kbhi, lblo, lbhi 534 double precision ssbb(lblo:lbhi,kblo:kbhi,jlo:jhi,ilo:ihi) 535 integer reclen 536 integer lwidth 537 integer iolab(lwidth,reclen) 538 double precision ioval(reclen) 539 integer iobuflen 540 double precision iobuf(*) 541 542 integer i, j, k, l 543 integer iir, issparse, lab1, lab2 544 integer recnum 545 double precision xx 546#include "bitops_decls.fh" 547#include "bitops_funcs.fh" 548 549 issparse = 1 550 recnum = 0 551 iir = 0 552 call ifill( lwidth*reclen, 0, iolab, 1 ) 553 call dfill( reclen, 0.d0, ioval, 1 ) 554 do i=ilo,ihi 555 do j=jlo,jhi 556 lab1 = ior(ishft(i,16),j) 557 do k=kblo,kbhi 558 do l=lblo,lbhi 559 xx = ssbb(l,k,j,i) 560 if (abs(xx).gt.1.d-12) then 561 iir = iir + 1 562 ioval(iir) = xx 563 lab2 = ior(ishft(k,16),l) 564 if (lwidth.eq.1) then 565 iolab(1,iir) = ior(ishft(lab1,32),lab2) 566 elseif (lwidth.eq.2) then 567 iolab(1,iir) = lab1 568 iolab(2,iir) = lab2 569 endif 570 if (iir.eq.reclen) then 571 call moints_iorec_flush(fd, fptr, tasknum, iir, 572 $ issparse, blkid, iobuflen, 573 $ iobuf ) 574 recnum = recnum + 1 575 call ifill( lwidth*reclen, 0, iolab, 1 ) 576 call dfill( reclen, 0.d0, ioval, 1 ) 577 iir = 0 578 endif 579 endif 580 enddo 581 enddo 582 enddo 583 enddo 584 if (iir.gt.0) then 585 call moints_iorec_flush( fd, fptr, tasknum, iir, issparse, 586 $ blkid, iobuflen, iobuf ) 587 recnum = recnum + 1 588 iir = 0 589 endif 590 591 return 592 end 593 594 595 596 subroutine moints_aodisk_iorec2sprs_old( ilo, ihi, jlo, jhi, 597 $ kblo, kbhi, lblo, lbhi, 598 $ ssbb, reclen, lwidth, iolab, 599 $ ioval ) 600 implicit none 601 integer ilo, ihi, jlo, jhi 602 integer kblo, kbhi, lblo, lbhi 603 double precision ssbb(lblo:lbhi,kblo:kbhi,jlo:jhi,ilo:ihi) 604 integer reclen 605 integer lwidth 606 integer iolab(lwidth,reclen) 607 double precision ioval(reclen) 608 integer q, i, j, k, l 609 integer i16mask 610 integer onbitmask 611 external onbitmask 612#include "bitops_decls.fh" 613#include "bitops_funcs.fh" 614 615 i16mask = onbitmask(16) 616 if (lwidth.eq.1) then 617 do q=1,reclen 618 i = iand(ishft(iolab(1,q),-48),i16mask) 619 j = iand(ishft(iolab(1,q),-32),i16mask) 620 k = iand(ishft(iolab(1,q),-16),i16mask) 621 l = iand(iolab(1,q),i16mask) 622 ssbb(l,k,j,i) = ioval(q) 623 enddo 624 else if (lwidth.eq.2) then 625 do q=1,reclen 626 i = iand(ishft(iolab(1,q),-16),i16mask) 627 j = iand(iolab(1,q),i16mask) 628 k = iand(ishft(iolab(2,q),-16),i16mask) 629 l = iand(iolab(2,q),i16mask) 630 ssbb(l,k,j,i) = ioval(q) 631 enddo 632 endif 633 return 634 end 635 636 637 638 639 640 641 subroutine moints_iorec_flush( fd, fptr, tasknum, reccnt, 642 $ issparse, blkid, buflen, iobuf ) 643 implicit none 644#include "mafdecls.fh" 645#include "eaf.fh" 646#include "mointsmoaodef.fh" 647 integer fd 648 double precision fptr 649 integer tasknum 650 integer reccnt 651 integer issparse 652 integer blkid 653 integer buflen 654 double precision iobuf(buflen) 655c 656 integer blkinfo(MOINTS_NBLKINFO) 657 integer stat, bufbytes 658c 659 call ifill( MOINTS_NBLKINFO, 0, blkinfo, 1 ) 660 blkinfo(1) = tasknum 661 blkinfo(2) = reccnt 662 blkinfo(3) = issparse 663 blkinfo(4) = blkid 664 call icopy( MOINTS_NBLKINFO, blkinfo, 1, iobuf, 1 ) 665 bufbytes = ma_sizeof(MT_DBL,buflen,MT_BYTE) 666 stat = eaf_write(fd, fptr, iobuf, bufbytes ) 667 fptr = fptr + bufbytes 668c 669c$$$ WRITE(6,772) blkid 670c$$$ 772 FORMAT('---->',I8) 671c 672 673 return 674 end 675 676 677 678 679 680 681 logical function moints_iorec_next( blkid, buidstr ) 682 implicit none 683#include "errquit.fh" 684#include "mafdecls.fh" 685#include "eaf.fh" 686#include "mointsmoaodef.fh" 687#include "cmointsmoao.fh" 688 integer blkid 689 character*8 buidstr 690 integer blkinfo(MOINTS_NBLKINFO) 691 integer stat 692 integer bufbytes 693 694 moints_iorec_next = .false. 695 if (moao_eof) return 696c 697c check if io buffer is already pre-read 698c 699 if (.not.(moao_touch)) then 700 if (blkid.eq.moao_blkid) then 701 moao_touch = .true. 702 moints_iorec_next = .true. 703 endif 704 return 705 endif 706c 707c otherwise read in a buffer 708c 709 bufbytes = ma_sizeof(MT_DBL,moao_buflen,MT_BYTE) 710 stat = eaf_read( moao_fd, moao_fptr, dbl_mb(moao_kbuf), 711 $ bufbytes ) 712 moao_fptr = moao_fptr + bufbytes 713 MOAO_NREC = MOAO_NREC + 1 714c 715c check for EOF and/or update info 716c 717 if (stat.eq.0) then 718 call icopy( MOINTS_NBLKINFO, dbl_mb(moao_kbuf), 1, blkinfo, 1 ) 719 moao_reclen = blkinfo(2) 720 moao_issparse = blkinfo(3) 721 moao_touch = (moao_tasknum.eq.blkinfo(1)).and. 722 $ (blkid.eq.blkinfo(4)) 723 if (moao_touch) then 724 moints_iorec_next = .true. 725 return 726 endif 727 moao_tasknum = blkinfo(1) 728 moao_blkid = blkinfo(4) 729c$$$ WRITE(6,772) moao_blkid 730c$$$ 772 FORMAT('<----',I8) 731 elseif (stat.gt.0) then 732 call errquit('moints_io disk io error',0, DISK_ERR) 733 endif 734c 735c reach here if tasknumber has changed or EOF 736c return FALSE 737c 738 moao_eof = (stat.lt.0) 739 return 740 end 741 742 743c 744c Create a unique identifier string for 745c 4 shell labels - 2 bytes per label 746c 747 subroutine moints_uid( str, i, j, k, l) 748 implicit none 749 character*8 str 750 integer i, j, k, l 751 character*2 ci, cj, ck, cl 752 integer*2 ii, jj, kk, ll 753 equivalence (ci,ii),(cj,jj),(ck,kk),(cl,ll) 754 755 ii = i 756 jj = j 757 kk = k 758 ll = l 759 str(1:2) = ci 760 str(3:4) = cj 761 str(5:6) = ck 762 str(7:8) = cl 763 764 return 765 end 766 767 768 769 770c 771c end OLD_AODISK section 772c 773#endif 774 775 776 777 778c 779c ============================================== 780c ============================================== 781c 782c 783c 784c New Version 785c 786c 787c 788c ============================================== 789c ============================================== 790c 791c 792 793 794 795 796 797#ifndef OLD_AODISK 798c 799c Alternate version of sparse packing 800c Labels/indices require integer of storage, 2*NCOL + NNZ/2, 801c cf. 2*NNZ for standard version. 802c 803c 804c IO Record structure: 805c The first integer of each IO record points to the 806c (double word) location of the first header. The space 807c between is assumed to the continuation of the 808c previous entry. ie. 809c 810c __________________________ 811c | | 812c | v 813c -------------------------------------------- 814c | | header.... 815c -------------------------------------------- 816c |<-- continuation -->| 817c prev record 818c 819c 820c The header structure is: 821c 822c hdr(1) magic I valid header 823c magic II invalid header, read next record (rare) 824c 825c hdr(2) >0 location of next header 826c -1 last header for this IO record 827c 828c 829c 830 831 subroutine moints_aodisk_sprs2iorec( fd, fptr, tasknum, blkid, 832 $ ilo, ihi, jlo, jhi, 833 $ kblo, kbhi, lblo, lbhi, 834 $ ssbb, recptr, reclen, iorec, 835 $ nrec ) 836 implicit none 837#include "errquit.fh" 838#include "mafdecls.fh" 839#include "mointsmoaodef.fh" 840 integer fd 841 double precision fptr 842 integer tasknum 843 integer blkid 844 integer ilo, ihi, jlo, jhi 845 integer kblo, kbhi, lblo, lbhi 846 double precision ssbb(lblo:lbhi,kblo:kbhi,jlo:jhi,ilo:ihi) 847 integer recptr 848 integer reclen 849 double precision iorec(reclen) 850 integer nrec 851c 852 integer i, j, ilen, jlen, llen, klen 853 integer swordlen, nir, nnz, dbilen, dijlen, nir1, nnz1 854 integer rp, hdrp, ijp, rleft, bblen, nssbb 855 integer irp, irrp, drp, hoff 856 integer ijidx, issparse 857 integer l_ij, k_ij 858 integer blkinfo(MOINTS_NBLKINFO) 859 integer iminus, magic_cookie, magic_cookie2 860 logical iscont, ojustfit, ohdrwrite, ojust1blk 861 integer labint 862 double precision sprstol 863 integer nblk,nexpand 864 double precision sparsesize, densesize 865 common/aodisk_stat/nblk,nexpand,sparsesize,densesize 866c 867 integer onbitmask 868 external onbitmask 869 integer dnonzero_cnt 870 external dnonzero_cnt 871c 872#include "bitops_decls.fh" 873 874*:: data statements must come after all declarations 875 876 data issparse/1/ 877 data iminus/-1/ 878 data sprstol/1.d-12/ 879 880#include "bitops_funcs.fh" 881c 882 blkinfo(1) = -1 883 blkinfo(3) = 0 884 blkinfo(4) = tasknum 885 blkinfo(5) = blkid 886 blkinfo(6) = issparse 887 blkinfo(7) = ilo 888 blkinfo(8) = ihi 889 blkinfo(9) = jlo 890 blkinfo(10) = jhi 891 blkinfo(11) = kblo 892 blkinfo(12) = kbhi 893 blkinfo(13) = lblo 894 blkinfo(14) = lbhi 895 ohdrwrite = .false. 896c 897c 898C WRITE(6,910) BLKINFO(4), BLKINFO(5), BLKINFO(7), 899C $ BLKINFO(8), BLKINFO(9), BLKINFO(10), 900C $ BLKINFO(11), BLKINFO(12), BLKINFO(13), 901C $ BLKINFO(14) 902C 910 FORMAT('Task#:',I5,5x,' ID =',I5,3x,4(i2,'-',i2,3x)) 903c 904 llen = lbhi - lblo + 1 905 klen = kbhi - kblo + 1 906 ilen = ihi - ilo + 1 907 jlen = jhi - jlo + 1 908 swordlen = ma_sizeof(MT_DBL,1,MT_BYTE) + 2 ! double word + 2 byte label = 10 909 magic_cookie = onbitmask(17) 910 magic_cookie2 = onbitmask(21) 911 labint = ma_sizeof(MT_INT,1,MT_BYTE)/2 ! labels per integer word, 2 bytes per label 912c 913c Minimium header length must be remaining on IO record 914c Otherwise flush IO record 915c 916 dbilen = ma_sizeof(MT_INT,MOINTS_NBLKINFO,MT_DBL) 917 dijlen = ma_sizeof(MT_INT,(3*ilen*jlen),MT_DBL) 918 rleft = reclen - recptr + 1 919 if ((dbilen+dijlen).gt.rleft) then 920 call moints_aodisk_flushiorec( fd, fptr, reclen, iorec ) 921 recptr = 2 922 call icopy( 1, recptr, 1, iorec(1), 1) 923 endif 924c 925c Header and record pointers 926c 927 ojust1blk = recptr.eq.2 928 hdrp = recptr 929 ijp = hdrp + dbilen 930 rp = hdrp + dbilen + dijlen 931 nssbb = 0 932c 933c Allocate temp space 934c 935 if (.not. ma_push_get(MT_INT, (ilen*jlen), 'ij ',l_ij, k_ij)) 936 $ call errquit('moints_aodisk_sprs2iorec: no memory',0, MA_ERR) 937 call ifill((ilen*jlen), 0, int_mb(k_ij), 1 ) 938c 939c 940c 941 iscont = .false. 942 do i=ilo,ihi 943 do j=jlo,jhi 944 945c 946c Count non-zeros -- compute space requirements 947c 948 rleft = reclen - rp + 1 949 nnz = dnonzero_cnt((llen*klen),sprstol,ssbb(lblo,kblo,j,i)) 950 nir = nnz/labint 951 if (mod(nnz,labint).ne.0) nir = nir + 1 952 bblen = ma_sizeof(MT_INT,(2*klen),MT_DBL) + 953 $ ma_sizeof(MT_INT,nir,MT_DBL) + nnz 954c 955c Compression statistics 956c 957 nblk = nblk + 1 958 if (bblen.gt.(klen*llen)) nexpand = nexpand + 1 959 sparsesize = sparsesize + bblen 960 densesize = densesize + llen*klen 961c 962c Flush IO record if insufficient space 963c 964 965 if (bblen.gt.rleft) then 966 if (.not.(ohdrwrite)) then 967 blkinfo(1) = magic_cookie 968 blkinfo(2) = -1 ! flag -- last entry in this record 969 call icopy( MOINTS_NBLKINFO, blkinfo, 1, iorec(hdrp), 1 ) 970 ohdrwrite = .true. 971 endif 972 call icopy( (ilen*jlen), int_mb(k_ij), 1, iorec(ijp), 1 ) 973 call moints_aodisk_flushiorec( fd, fptr, reclen, iorec ) 974c 975c Reset pointers and arrays 976c 977 call icopy( 1, iminus, 1, iorec(1), 1) 978 call ifill((ilen*jlen), 0, int_mb(k_ij), 1 ) 979 ijp = 2 980 rp = ijp + dijlen 981 nssbb = 0 982 iscont = .true. 983 ojust1blk = .true. 984 endif 985c 986c Pack sparse structure into IO record 987c 988 ijidx = (i-ilo)*jlen + j - jlo 989 int_mb(k_ij+ijidx) = ior(ishft(rp,16),nnz) 990 nssbb = nssbb + 1 991 irp = rp 992 irrp = irp + ma_sizeof(MT_INT,(2*klen),MT_DBL) 993 drp = irrp + ma_sizeof(MT_INT,nir,MT_DBL) 994 call moints_sparse2d_pack( llen, klen, ssbb(lblo,kblo,j,i), 995 $ sprstol, labint, iorec(irp), nir1, 996 $ iorec(irrp), iorec(drp), nnz1 ) 997 rp = drp + nnz1 998 if ((nnz1.ne.nnz).or.(nir1.ne.nir)) 999 $ call errquit('moints_aodisk_sprs2iorec: internal error 2',0, 1000 & DISK_ERR) 1001 enddo 1002 1003 enddo 1004 if (rp.gt.(reclen+1)) 1005 $ call errquit('moints_aodisk_sprs2iorec: internal error 1',0, 1006 & DISK_ERR) 1007c 1008c Bookkeeping for some odd cases 1009c 1010 ojustfit = (rp.eq.(reclen+1)) ! special case 1011 if (ojustfit) then 1012 rp = 0 1013 if (ojust1blk) then 1014 hoff = 0 1015 if (.not.(ohdrwrite)) hoff = hdrp 1016 call icopy( 1, hoff, 1, iorec(1), 1 ) 1017 endif 1018 endif 1019c 1020c Copy header info & map to IO record before returning 1021c 1022 if (.not.(ohdrwrite)) then 1023 blkinfo(1) = magic_cookie 1024 blkinfo(2) = rp 1025 call icopy( MOINTS_NBLKINFO, blkinfo, 1, iorec(hdrp), 1 ) 1026 endif 1027 call icopy( (ilen*jlen), int_mb(k_ij), 1, iorec(ijp), 1 ) 1028c 1029c Special case -- flush IO record 1030c 1031 if (ojustfit) then 1032 call moints_aodisk_flushiorec( fd, fptr, reclen, iorec) 1033 rp = 2 1034 call icopy( 1, rp, 1, iorec(1), 1) 1035 endif 1036c 1037c Put magic cookie II at start of next header. 1038c Usually overwritten by magic cookie next time around 1039c (unless we skip to next record) 1040c 1041 if (rp.le.reclen) then 1042 call icopy(1, magic_cookie2, 1, iorec(rp), 1 ) 1043 endif 1044c 1045c If SSBB block spans multiple records must encode offset 1046c of end of block 1047c 1048 if (iscont) then 1049 call icopy( 1, rp, 1, iorec(1), 1 ) 1050 endif 1051 recptr = rp 1052c 1053c Clean up 1054c 1055 if (.not. ma_pop_stack(l_ij)) 1056 $ call errquit('moints_aodisk_sprs2iorec: failed to pop', l_ij, 1057 & MA_ERR) 1058 1059 return 1060 end 1061 1062 1063 1064 1065 1066 1067c 1068c Read and interpret AO file contents - one block 1069c at a time -- structure as above 1070c 1071c 1072 subroutine moints_aodisk_iorec2sprs( fd, fptr, tasknum, 1073 $ ilo, ihi, jlo, jhi, 1074 $ kblo, kbhi, lblo, lbhi, 1075 $ ssbb, hdrp, reclen, iorec, 1076 $ nrec ) 1077 implicit none 1078#include "errquit.fh" 1079#include "mafdecls.fh" 1080#include "eaf.fh" 1081#include "mointsmoaodef.fh" 1082 integer fd 1083 double precision fptr 1084 integer tasknum 1085 integer ilo, ihi, jlo, jhi 1086 integer kblo, kbhi, lblo, lbhi 1087 double precision ssbb(lblo:lbhi,kblo:kbhi,jlo:jhi,ilo:ihi) 1088 integer hdrp 1089 integer reclen 1090 double precision iorec(reclen) 1091 integer nrec 1092c 1093 integer cookie, cookie2 1094 integer dbilen, ilen, jlen, ijp 1095 integer nexthdrp, stat, bufbytes, magic 1096 INTEGER IJ 1097 integer l_ij, k_ij 1098 integer blkinfo(MOINTS_NBLKINFO) 1099 logical ocont 1100 integer onbitmask 1101 external onbitmask 1102c 1103c 1104c 1105 cookie = onbitmask(17) 1106 cookie2 = onbitmask(21) 1107 bufbytes = ma_sizeof(MT_DBL,reclen,MT_BYTE) 1108 ilen = ihi - ilo + 1 1109 jlen = jhi - jlo + 1 1110 dbilen = ma_sizeof(MT_INT,MOINTS_NBLKINFO,MT_DBL) 1111 ijp = hdrp + dbilen 1112 if (.not. ma_push_get(MT_INT, (ilen*jlen), 'ij ',l_ij, k_ij)) 1113 $ call errquit('moints_aodisk_iorec2sprs: no memory',0, 1114 & MA_ERR) 1115c 1116c Sanity check on current pointer and IO record 1117c 1118 call icopy( MOINTS_NBLKINFO, iorec(hdrp), 1, blkinfo, 1) 1119 if (blkinfo(1).ne.cookie) 1120 $ call errquit('moints_aodisk_iorec2sprs: AOfile corrupt I?',0, 1121 & DISK_ERR) 1122 1123c$$$ WRITE(6,910) BLKINFO(4), BLKINFO(5), BLKINFO(7), 1124c$$$ $ BLKINFO(8), BLKINFO(9), BLKINFO(10), 1125c$$$ $ BLKINFO(11), BLKINFO(12), BLKINFO(13), 1126c$$$ $ BLKINFO(14) 1127c$$$ 910 FORMAT('Task#:',I5,5x,' ID =',I5,3x,4(i2,'-',i2,3x)) 1128c 1129c Check consistency between args and blockinfo 1130c 1131 if ((ilo.ne.blkinfo(7)).or.(ihi.ne.blkinfo(8)).or. 1132 $ (jlo.ne.blkinfo(9)).or.(jhi.ne.blkinfo(10)).or. 1133 $ (kblo.ne.blkinfo(11)).or.(kbhi.ne.blkinfo(12)).or. 1134 $ (lblo.ne.blkinfo(13)).or.(lbhi.ne.blkinfo(14))) then 1135 write(6,623) nrec,(blkinfo(ij),ij=7,14) 1136 623 format(' Rec#',i6,3x,' B ',8i6) 1137 write(6,622) nrec,ilo,ihi,jlo,jhi,kblo,kbhi,lblo,lbhi 1138 622 format(' Rec#',i6,3x,' E ',8i6) 1139 call errquit('moints_aodisk_sprs2iorec: blk label mismatch',0, 1140 & INPUT_ERR) 1141 endif 1142c 1143c Unpack data structure from IO record(s) 1144c -- possibly spanning multiple records 1145c 1146 ocont = .true. 1147 nexthdrp = blkinfo(2) 1148 1149 do while (ocont) 1150 1151 call icopy( (ilen*jlen), iorec(ijp), 1, int_mb(k_ij), 1) 1152 1153c$$$ CALL MOINTS_AODISK_IJPRINT( ILEN, JLEN, INT_MB(K_IJ) ) 1154 1155 call moints_aodisk_sprs2dense_a( ilo, ihi, jlo, jhi, 1156 $ kblo, kbhi, lblo, lbhi, 1157 $ int_mb(k_ij), reclen, iorec, 1158 $ ssbb ) 1159 ocont = nexthdrp.eq.-1 1160 if (nexthdrp.le.0) then 1161 stat = eaf_read(fd, fptr, iorec, bufbytes ) 1162 nrec = nrec + 1 1163 if (stat.ne.0) goto 150 1164 fptr = fptr + bufbytes 1165 call icopy( 1, iorec(1), 1, nexthdrp, 1 ) 1166 ijp = 2 1167 endif 1168 enddo 1169c 1170c Check next block -- skip ahead if invalid 1171c 1172 101 call icopy( 1, iorec(nexthdrp), 1, magic, 1) 1173 if (magic.eq.cookie2) then 1174 stat = eaf_read(fd, fptr, iorec, bufbytes ) 1175 nrec = nrec + 1 1176 if (stat.ne.0) goto 150 1177 fptr = fptr + bufbytes 1178 nexthdrp = 2 1179 goto 101 1180 endif 1181 if (magic.ne.cookie) 1182 $ call errquit('moints_aodisk_iorec2sprs: AOfile corruptII?',0, 1183 & DISK_ERR) 1184c 1185c Prefetch tasknumber from header 1186c 1187 call icopy( 4, iorec(nexthdrp), 1, blkinfo, 1 ) 1188 tasknum = blkinfo(4) 1189 hdrp = nexthdrp 1190c 1191c Cleanup 1192c 1193 150 continue 1194 if (.not. ma_pop_stack(l_ij)) 1195 $ call errquit('moints_aodisk_iorec2sprs: failed to pop', l_ij, 1196 & MA_ERR) 1197 return 1198 end 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 subroutine moints_aodisk_flushiorec( fd, fptr, n, r ) 1210 implicit none 1211#include "mafdecls.fh" 1212#include "eaf.fh" 1213#include "cmointsmoao.fh" 1214 integer fd 1215 double precision fptr 1216 integer n 1217 double precision r(n) 1218 integer bufbytes 1219 integer stat 1220 1221 moao_nrec = moao_nrec + 1 1222 bufbytes = ma_sizeof(MT_DBL,n,MT_BYTE) 1223 stat = eaf_write(fd, fptr, r, bufbytes ) 1224 fptr = fptr + bufbytes 1225 return 1226 end 1227 1228 1229 1230 1231 1232 subroutine moints_aodisk_sprs2dense_a( ilo, ihi, jlo, jhi, 1233 $ kblo, kbhi, lblo, lbhi, 1234 $ ijmap, reclen, iorec, ssbb) 1235 implicit none 1236#include "errquit.fh" 1237#include "mafdecls.fh" 1238 integer ilo, ihi, jlo, jhi 1239 integer kblo, kbhi, lblo, lbhi 1240 integer ijmap(jlo:jhi,ilo:ihi) 1241 integer reclen 1242 double precision iorec(*) 1243 double precision ssbb(lblo:lbhi,kblo:kbhi,jlo:jhi,ilo:ihi) 1244c 1245 integer labint 1246 integer i, j, ijp, nnz, bit16, ilen, jlen, klen, llen 1247 integer maxiirlen, iirlen, nir, dp 1248 integer l_iv, k_iv, k_irv 1249 integer onbitmask 1250 external onbitmask 1251#include "bitops_decls.fh" 1252#include "bitops_funcs.fh" 1253 1254 ilen = ihi - ilo + 1 1255 jlen = jhi - jlo + 1 1256 klen = kbhi - kblo + 1 1257 llen = lbhi - lblo + 1 1258 bit16 = onbitmask(16) 1259 maxiirlen = 2*klen + (klen*llen)/2 + 1 1260 if (.not. ma_push_get(MT_INT, maxiirlen, 'i io pack ',l_iv, k_iv)) 1261 $ call errquit('moints_aodisk_sprs2dense_a: no memory',0, MA_ERR) 1262 k_irv = k_iv + 2*klen 1263 labint = ma_sizeof(MT_INT,1,MT_BYTE)/2 1264 1265 do i=ilo,ihi 1266 do j=jlo,jhi 1267 ijp = iand(ishft(ijmap(j,i),-16),bit16) 1268 nnz = iand(ijmap(j,i),bit16) 1269 if (nnz.gt.(klen*llen)) then 1270 call moints_aodisk_dumpiorec( reclen, iorec ) 1271 call errquit('moints_sprs2dense_a: internal error 3',0, 1272 & INT_ERR) 1273 endif 1274 if (ijp.gt.0) then 1275 nir = nnz/labint 1276 if (mod(nnz,labint).ne.0) nir = nir + 1 1277 iirlen = 2*klen + nir 1278 if (iirlen.gt.maxiirlen) then 1279 stop 1190 1280 endif 1281 dp = ijp + ma_sizeof(MT_INT,iirlen,MT_DBL) 1282 call icopy( iirlen, iorec(ijp), 1, int_mb(k_iv), 1) 1283 call moints_sparse2d_unpack( llen, klen, int_mb(k_iv), 1284 $ nnz, int_mb(k_irv), iorec(dp), 1285 $ ssbb(lblo,kblo,j,i) ) 1286 endif 1287 enddo 1288 enddo 1289 if (.not. ma_pop_stack(l_iv)) 1290 $ call errquit('moints_aodisk_sprs2iorec: failed to pop', l_iv, 1291 & MA_ERR) 1292 1293 end 1294 1295c 1296c Define this macro assume everyone supports integer*2 1297c Undefine if compiler complains 1298c 1299#define HAVE_FORTRAN_INTEGER2 1300 1301c 1302c This routine packs a dense 2d matrix (m x n) 1303c into a standard sparse structure, 1304c column pointer with non-zero rows, row-indices 1305c 1306c 1307c COLPLO[]: 1 2 3 4 5 .... n 1308c | | 1309c | ----------- 1310c | | 1311c V V 1312c V[]: 1 .... 9 ... 1313c IR[]: 2 4 5 8... 2 4 9 1314c 1315c 1316c 1317#ifdef HAVE_FORTRAN_INTEGER2 1318 subroutine moints_sparse2d_pack( m, n, x, tol, labint, 1319 $ colp, nir, ir, v, nnz ) 1320 implicit none 1321 integer m ! [input] rows 1322 integer n ! [input] columns 1323 double precision x(m,n) ! [input] dense 2D 1324 double precision tol ! [input] tolerance 1325 integer labint ! [input] 2byte labels per integer word 1326 integer colp(2,n) ! [output] column ptr hi,lo 1327 integer nir ! [output] row index length 1328 integer*2 ir(*) ! [output] row index 16-bit packed 1329 double precision v(*) ! [output] packed values 1330 integer nnz ! [output] number non-zeroes 1331c 1332 integer s, t 1333 double precision xx 1334 1335 1336 nnz = 0 1337 do s=1,n 1338 colp(1,s) = nnz + 1 1339 do t=1,m 1340 xx = x(t,s) 1341 if (abs(xx).ge.tol) then 1342 nnz = nnz + 1 1343 v(nnz) = xx 1344 ir(nnz) = t 1345 endif 1346 enddo 1347 if (colp(1,s).gt.nnz) colp(1,s) = 0 1348 colp(2,s) = nnz 1349 enddo 1350 nir = nnz/labint 1351 if (mod(nnz,labint).ne.0) nir = nir + 1 1352 1353 return 1354 end 1355 1356 1357 1358 1359 1360 subroutine moints_sparse2d_unpack( m, n, colp, nnz, ir, v, x ) 1361 implicit none 1362 integer m ! [input] rows 1363 integer n ! [input] columns 1364 integer colp(2,n) ! [input] column ptr hi,lo 1365 integer nnz ! [input] number non-zeroes 1366 double precision v(nnz) ! [input] packed values 1367 double precision x(m,n) ! [output] dense 2D 1368 integer*2 ir(*) ! [input] row index 16-bit packed 1369 integer s, rp 1370 1371 do s=1,n 1372 if (colp(1,s).gt.0) then 1373 do rp=colp(1,s),colp(2,s) 1374 x(ir(rp),s) = v(rp) 1375 enddo 1376 endif 1377 enddo 1378 return 1379 end 1380 1381 1382 1383c 1384c ---------------------------------------------- 1385c Alternative code without integer*2 data type 1386c ---------------------------------------------- 1387c 1388#else 1389 subroutine moints_sparse2d_pack( m, n, x, tol, labint, colp, 1390 $ nir, ir, v, nnz ) 1391 implicit none 1392 integer m ! [input] rows 1393 integer n ! [input] columns 1394 double precision x(m,n) ! [input] dense 2D 1395 double precision tol ! [input] tolerance 1396 integer labint ! [input] 2byte labels per integer word 1397 integer colp(2,n) ! [output] column ptr hi,lo 1398 integer nir ! [output] row index length 1399 integer ir(*) ! [output] row index 16-bit packed 1400 double precision v(*) ! [output] packed values 1401 integer nnz ! [output] number non-zeroes 1402c 1403 integer s, t, bp, pack 1404 integer ibit16 1405 double precision xx 1406#include "bitops_decls.fh" 1407#include "bitops_funcs.fh" 1408 integer onbitmask 1409 external onbitmask 1410 1411 ibit16 = onbitmask(16) 1412 nnz = 0 1413 nir = 0 1414 pack = 0 1415 bp = 0 1416 do s=1,n 1417 colp(1,s) = nnz + 1 1418 do t=1,m 1419 xx = x(t,s) 1420 if (abs(xx).gt.tol) then 1421 nnz = nnz + 1 1422 v(nnz) = xx 1423 bp = mod((labint-mod(nnz,labint)),labint)*16 1424C bp = mod((2-mod(nnz,2)),2)*16 1425 pack = ior(pack,ishft(t,bp)) 1426 if (bp.eq.0) then 1427 nir = nir + 1 1428 ir(nir) = pack 1429 pack = 0 1430 if (nnz.ne.(labint*nir)) stop 3331 1431 endif 1432 endif 1433 enddo 1434 if (colp(1,s).gt.nnz) colp(1,s) = 0 1435 colp(2,s) = nnz 1436 enddo 1437 if (bp.ne.0) then 1438 nir = nir + 1 1439 ir(nir) = pack 1440 endif 1441 return 1442 end 1443 1444 1445 1446 1447 1448 subroutine moints_sparse2d_unpack( m, n, colp, nnz, ir, v, x ) 1449 implicit none 1450#include "mafdecls.fh" 1451 integer m ! [input] rows 1452 integer n ! [input] columns 1453 integer colp(2,n) ! [input] column ptr hi,lo 1454 integer nnz ! [input] number non-zeroes 1455 integer ir(*) ! [input] row index 16-bit packed 1456 double precision v(nnz) ! [input] packed values 1457 double precision x(m,n) ! [output] dense 2D 1458 integer s, t, rp, rrp, parity, nshft 1459 integer bit16, labint 1460#include "bitops_decls.fh" 1461#include "bitops_funcs.fh" 1462 integer onbitmask 1463 external onbitmask 1464 1465 labint = ma_sizeof(MT_INT,1,MT_BYTE)/2 1466 bit16 = onbitmask(16) 1467 do s=1,n 1468 if (colp(1,s).gt.0) then 1469 do rp=colp(1,s),colp(2,s) 1470 parity = mod(rp,labint) 1471 rrp = rp/labint 1472 if (parity.ne.0) rrp = rrp + 1 1473 nshft = 16*(labint - parity) 1474 if (parity.eq.0) nshft = 0 1475 t = iand(ishft(ir(rrp),-nshft),bit16) 1476 x(t,s) = v(rp) 1477 enddo 1478 endif 1479 enddo 1480 return 1481 end 1482 1483#endif 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 integer function dnonzero_cnt( n, tol, v ) 1496 implicit none 1497 integer n 1498 double precision tol 1499 double precision v(n) 1500 integer i, nnz 1501 1502 nnz = 0 1503 do i=1,n 1504 if (abs(v(i)).ge.tol) nnz = nnz + 1 1505 enddo 1506 dnonzero_cnt = nnz 1507 return 1508 end 1509 1510 1511 1512 1513 1514c 1515c ========================================================= 1516c 1517c Debugging and diagnostic codes 1518c 1519c ========================================================= 1520c 1521 1522 1523 1524 1525c 1526c Debugging code to check file consistency 1527c 1528 subroutine moints_aodisk_checkfile( fd, reclen, iorec ) 1529 implicit none 1530#include "mafdecls.fh" 1531#include "eaf.fh" 1532#include "mointsmoaodef.fh" 1533 integer fd 1534 integer reclen 1535 double precision iorec(reclen) 1536 integer stat, recbytes, rp, hdrp, ijp 1537 integer ilen, jlen 1538 integer dbilen 1539 double precision fptr 1540 integer blkinfo(MOINTS_NBLKINFO) 1541 integer magic_cookie, magic_cookie2 1542 integer ijmap(1000) 1543 integer onbitmask 1544 external onbitmask 1545 1546 magic_cookie = onbitmask(17) 1547 magic_cookie2 = onbitmask(21) 1548 dbilen = ma_sizeof(MT_INT,MOINTS_NBLKINFO,MT_DBL) 1549 stat = 0 1550 recbytes = ma_sizeof(MT_DBL,reclen,MT_BYTE) 1551 fptr = 1.d0 1552 rp = 1 1553 hdrp = 2 1554 blkinfo(2) = 0 1555 ilen = 0 1556 jlen = 0 1557 1558 100 stat = eaf_read(fd, fptr, iorec, recbytes ) 1559 if (stat.eq.0) then 1560 fptr = fptr + recbytes 1561 if ((blkinfo(2).eq.-1).and.(hdrp.ne.0)) then 1562 call icopy((ilen*jlen), iorec(2), 1, ijmap, 1 ) 1563C call moints_aodisk_ijprint( ilen, jlen, ijmap ) 1564 endif 1565 1566 call icopy(1, iorec(1), 1, hdrp, 1 ) 1567 if (hdrp.le.0) goto 100 1568 101 call icopy(MOINTS_NBLKINFO, iorec(hdrp), 1, blkinfo, 1 ) 1569 ijp = hdrp + dbilen 1570 if (blkinfo(1).eq.magic_cookie2) then 1571c$$$ write(6,912) 1572c$$$ 912 format('header invalid *') 1573 goto 100 1574 endif 1575 if (blkinfo(1).ne.magic_cookie) stop 6669 1576 1577 1578 ilen = blkinfo(8) - blkinfo(7) + 1 1579 jlen = blkinfo(10) - blkinfo(9) + 1 1580 call icopy((ilen*jlen), iorec(ijp), 1, ijmap, 1 ) 1581c$$$ call moints_aodisk_ijprint( ilen, jlen, ijmap ) 1582 1583 if (blkinfo(2).gt.0) then 1584 hdrp = blkinfo(2) 1585 goto 101 1586 endif 1587 goto 100 1588 endif 1589 1590 return 1591 end 1592 1593 1594 1595 1596 subroutine moints_aodisk_dumpiorec( n, r ) 1597 implicit none 1598#include "mafdecls.fh" 1599#include "mointsmoaodef.fh" 1600 integer n 1601 double precision r(n) 1602 1603 integer hdrp, ijp, i, ilen, jlen 1604 integer bi(MOINTS_NBLKINFO) 1605 1606 call icopy(1, r(1), 1, hdrp, 1) 1607 do while (hdrp.gt.0) 1608 call icopy(MOINTS_NBLKINFO, r(hdrp), 1, bi, 1) 1609 write(6,901) (bi(i),i=1,MOINTS_NBLKINFO) 1610 901 format(16i8) 1611 ilen = bi(8) - bi(7) + 1 1612 jlen = bi(10) - bi(9) + 1 1613 ijp = hdrp + ma_sizeof(MT_INT,MOINTS_NBLKINFO,MT_DBL) 1614 call moints_aodisk_ijprint( ilen, jlen, r(ijp) ) 1615 1616 hdrp = bi(2) 1617 enddo 1618 return 1619 end 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 logical function moints_aodisk_verifyiorec(ilen, jlen, klen, 1631 $ llen, n, r ) 1632 implicit none 1633#include "mafdecls.fh" 1634#include "mointsmoaodef.fh" 1635 integer ilen, jlen, klen, llen 1636 integer n 1637 double precision r(n) 1638 1639 integer hdrp, cookie 1640 integer bi(MOINTS_NBLKINFO) 1641 integer itmp(1000) 1642 integer irec, ijp 1643 logical ao_ijverify 1644 external ao_ijverify 1645 integer onbitmask 1646 external onbitmask 1647 1648 moints_aodisk_verifyiorec = .true. 1649 cookie = onbitmask(17) 1650 irec = 1 1651 ijp = 2 1652 call icopy(1, r(1), 1, hdrp, 1) 1653 do while ((irec.eq.1).or.(hdrp.gt.0)) 1654 if ((hdrp.eq.2).or.(irec.gt.1)) then 1655 call icopy(MOINTS_NBLKINFO, r(hdrp), 1, bi, 1) 1656 if (bi(1).ne.cookie) return 1657 ijp = hdrp + ma_sizeof(MT_INT,MOINTS_NBLKINFO,MT_DBL) 1658 endif 1659 call icopy( (ilen*jlen), r(ijp), 1, itmp, 1) 1660 if (.not.(ao_ijverify( ilen, jlen, klen, llen, 1661 $ itmp, n, r))) then 1662 write(6,881) 1663 881 format('Failed to verify') 1664 moints_aodisk_verifyiorec = .false. 1665 return 1666 endif 1667 if ((irec.gt.1).or.(hdrp.eq.2)) hdrp = bi(2) 1668 irec = irec + 1 1669 enddo 1670 return 1671 end 1672 1673 1674 1675 1676 1677 1678 subroutine moints_aodisk_ijprint( jlen, ilen, ijmap ) 1679 implicit none 1680 integer ilen, jlen 1681 integer*2 ijmap(2,jlen,ilen) 1682 integer i, j 1683 1684 do j=1,jlen 1685 write(6,771) (ijmap(1,j,i),i=1,ilen) 1686 771 format(16i5) 1687 enddo 1688 return 1689 end 1690 1691 1692 1693 1694 1695 1696 logical function ao_ijverify( ilen, jlen, klen, llen, 1697 $ ijmap, n, r ) 1698 implicit none 1699 integer ilen, jlen, klen, llen 1700 integer*2 ijmap(2,jlen,ilen) 1701 integer n 1702 double precision r(n) 1703c 1704 integer i, j 1705 integer ijp, nnz, nir, iirlen 1706 integer itmp(1000) 1707 logical sprs_verify 1708 external sprs_verify 1709 1710 ao_ijverify = .true. 1711 do i=1,ilen 1712 do j=1,jlen 1713 if (ijmap(1,j,i).ne.0) then 1714 ijp = ijmap(1,j,i) 1715 nnz = ijmap(2,j,i) 1716 nir = nnz/2 + mod(nnz,2) 1717 iirlen = 2*klen + nir 1718 call icopy( iirlen, r(ijp), 1, itmp, 1) 1719 if (.not.(sprs_verify(llen,klen,itmp(1), 1720 $ itmp(2*klen+1)))) then 1721 print*,' --> CHECKING AT ',ijp 1722 write(6,771) i,j 1723 771 format('Failed to verify: I =',i5,' J =',I5) 1724 write(6,881) 1725 881 format(//,'IJ MAP dump') 1726 call moints_aodisk_ijprint(jlen,ilen,ijmap) 1727 ao_ijverify = .false. 1728 return 1729 endif 1730 endif 1731 enddo 1732 enddo 1733 return 1734 1735 end 1736 1737 1738 1739 1740 1741 1742 1743 subroutine sprs_print( m, n, colp, ir ) 1744 implicit none 1745 integer m ! [input] rows 1746 integer n ! [input] columns 1747 integer colp(2,n) ! [input] column ptr hi,lo 1748 integer*2 ir(*) 1749 integer s, rp 1750 1751 write(6,711) (colp(1,s),colp(2,s),s=1,n) 1752 711 format(//,16(i4,2x,i4,5x)) 1753 do s=1,n 1754 if (colp(1,s).gt.0) then 1755 write(6,772) s,(ir(rp),rp=colp(1,s),colp(2,s)) 1756 772 format(i5,5x,16i4) 1757 endif 1758 enddo 1759 return 1760 end 1761 1762 1763 logical function sprs_verify( m, n, colp, ir ) 1764 implicit none 1765 integer m ! [input] rows 1766 integer n ! [input] columns 1767 integer colp(2,n) ! [input] column ptr hi,lo 1768 integer*2 ir(*) 1769 integer s, t 1770 1771 sprs_verify = .true. 1772 do s=1,n 1773 if (colp(1,s).gt.0) then 1774 do t=colp(1,s),colp(2,s) 1775 if ((ir(t).le.0).or.(ir(t).gt.m)) then 1776 write(6,811) s,ir(t) 1777 811 format('Failed to verify: S=',i4,' T=',i4) 1778 call sprs_print( m, n, colp, ir ) 1779 sprs_verify = .false. 1780 return 1781 endif 1782 enddo 1783 endif 1784 enddo 1785 return 1786 end 1787 1788c 1789c End !OLD_AODISK section 1790c 1791#endif 1792