c$Id:$ c-----[--.----+----.----+----.-----------------------------------------] c File contains following USER modules: c (1) Subroutine umacr1: Command routine to call solver (USOL) c (2) Subroutine uasble: Assembly routine (called by DASBLE) c (3) Subroutine usolve: Solution driver (called by PRESOL & c PSOLVE) c Must have the SGI solution package which includes routines: c PSLDLT_ORDERING, PSLDLT_PREPROCESS, PSLDLT_FACTOR, and c PSDLT_SOLVE. c-----[--.----+----.----+----.-----------------------------------------] subroutine umacr1(lct,ctl,prt) c * * F E A P * * A Finite Element Analysis Program c.... Copyright (c) 1984-2009: Regents of the University of California c All Rights Reserved c-----[--.----+----.----+----.-----------------------------------------] c Purpose: Controls solution by SGI sparse solver. Also sets c storage for compressed arrays. c Command: USOL,,1 c Inputs: c lct - Command character parameters c ctl(3) - Command numerical parameters c prt - Flag, output if true c Outputs: c none - Outputs stored in blank common and pointers. c-----[--.----+----.----+----.-----------------------------------------] implicit none include 'allotd.h' include 'cdata.h' include 'comblk.h' include 'compac.h' include 'compas.h' include 'debugs.h' include 'iofile.h' include 'ndata.h' include 'part0.h' include 'pointer.h' include 'psize.h' include 'sdata.h' include 'setups.h' include 'ssolve.h' include 'umac1.h' real*8 opsa common /uslv/ opsa logical setvar,palloc, pcomp,prt character lct*15, tname*5 integer isw,j,kp,nonz real*4 tary(2), etime , tt real*8 ctl,ops save c Set command word if(pcomp(uct,'mac1',4)) then uct = 'usol' return endif if(pcomp(lct,'off',3)) then solver = .true. if(iow.lt.0) then write(*,*) ' System solver specified' endif else solver = .false. if(iow.lt.0) then write(*,*) ' User solver specified' endif diagin = .true. kbycol = .true. ! Store by columns c Compute sparse storage for non-zero matrix setvar = palloc(111,'TEMP1', numnp*ndf, 1) call elcnt(numel,nen,nen1,mr(np(31)),mr(np(33)),mr(np(111)),1) call sumcnt(mr(np(111)),numnp*ndf,kp) setvar = palloc(112,'TEMP2', kp, 1) call pelcon(numel,nen,nen1,mr(np(33)),mr(np(31)),mr(np(111)), & mr(np(112)),kp,1) setvar = palloc( 93,'OINC ', 2*neq+1 , 1) if(np(94).ne.0) then setvar = palloc( 94,'OINO ',0, 1) endif setvar = palloc( 94,'OINO ', 1, 1) call compro(numnp,nen,nen1,ndf,mr(np(33)),mr(np(31)), & mr(np(111)),mr(np(112)),mr(np(94)),mr(np(93)), & kp,maxm-mmax,kbycol,diagin,all) setvar = palloc( 94,'OINO ', kp, 1) c Delete temporary arrays setvar = palloc(112,'TEMP2', 0, 1) setvar = palloc(111,'TEMP1', 0, 1) c Set Harwell-Boeing pointer array into AD0 setvar = palloc( 67,'AD0 ', neq+1, 1) mr(np(67)) = 1 do j = 1,neq mr(np(67) + j) = mr(np(93)+j-1) + 1 end do ! j c Set storage for sparse stiffness array write(tname,'(4hTANG,i1)') npart setvar = palloc(npart,tname, kp, 2) call pzero(hr(np(npart)),kp) na = np(npart) nnr = kp nau = na + neq nal = nau numcels= 0 compfl = .true. c Sort row entries into increasing order call sortjc(mr(np(93)),mr(np(94)),neq) c Set 'p' and 'ip' to equation numbers setvar = palloc( 47, 'PNTER',neq, 1) setvar = palloc( 48, 'INVPT',neq, 1) do j = 1,neq mr(np(47)+j-1) = j mr(np(48)+j-1) = j end do ! j c Output solution properties tt = etime(tary) write(iow,2000) kp,tary if(ior.lt.0) then write(*,2000) kp,tary end if c Preprocess matrix to determine fill call PSLDLT_ORDERING(1,2) call PSLDLT_PREPROCESS(1, neq, mr(np(67)), mr(np(94)), & nonz, ops) if(ior.lt.0) then write(*,*) ' Number of non-zero entries =', nonz write(*,*) ' Number of factorization ops =', ops endif opsa = ops endif c Format 2000 format(10x,'Compressed Storage =',i9,20x,'t=',0p,2f9.2) end subroutine uasble(s,p,ld,ns,afl,bfl) c * * F E A P * * A Finite Element Analysis Program c.... Copyright (c) 1984-2009: Regents of the University of California c All Rights Reserved c-----[--+---------+---------+---------+---------+---------+---------+-] c Purpose: User assembly routine for SGI sparse solver c Inputs: c s(ns,ns) - element matrix c p(ns) - element vector c ld(ns) - local/global active equation numbers c ns - size of arrays c afl - Assemble s(ns,ns) into global storage c bfl - Assemble p(ns) into global storage c Outputs: c Assembled arrays in blank common. Address passed by pointer. c-----[--+---------+---------+---------+---------+---------+---------+-] implicit none include 'compac.h' include 'ndata.h' include 'pointer.h' include 'comblk.h' logical afl,bfl integer ns, nrhs, i, ld(ns) real*8 s(ns,ns),p(ns) save c Assemble compressed column with diagonals in if(afl) then call cassem(hr(na),hr(nau),hr(nal),s,mr(np(94)),mr(np(93)), & mr(np(48)),ld,ns,.false.,kbycol,diagin, all) endif c Assemble a vector if(bfl) then nrhs = np(26) - 1 do i = 1,ns if(ld(i).gt.0) hr(nrhs+ld(i)) = hr(nrhs+ld(i)) + p(i) end do endif end subroutine usolve(flags) c * * F E A P * * A Finite Element Analysis Program c.... Copyright (c) 1984-2009: Regents of the University of California c All Rights Reserved c-----[--+---------+---------+---------+---------+---------+---------+-] c Purpose: User solver interface c Inputs: c flags(1) - Allocation and/or initialization phase c flags(2) - Perform factorization for direct solutions c flags(3) - Coefficient array unsymmetric c flags(4) - Solve equations c flags(5) - .false. on entry c Outputs: c flags(5) - True if error occurs c-----[--+---------+---------+---------+---------+---------+---------+-] implicit none include 'cdata.h' include 'compas.h' include 'endata.h' include 'part0.h' include 'pointer.h' include 'comblk.h' real*8 ops common /uslv/ ops logical flags(*), setvar,palloc integer n real*4 etime, tt,tary(2),telps save c Initialize coefficient tangent array if(flags(1)) then call pzero(hr(np(npart)), nnr ) endif c Factor coefficient matrix: SGI sparse solver if(flags(2)) then tt = etime(tary) telps = tary(1) call psldlt_factor(1,neq,mr(np(67)),mr(np(94)),hr(np(npart))) tt = etime(tary) telps = tary(1) - telps write(*,*) 'USOLVE: Mflops = ', ops*1.d-6/dble(telps) endif c Solve equations: SGI sparse solver if(flags(4)) then setvar = palloc(111,'TEMP1',neq, 2) call psldlt_solve(1, hr(np(111)), hr(np(26))) write(*,*) ' Solve Equations in USOLVE ' aengy = 0.d0 do n = 0,neq-1 aengy = aengy + hr(np(111)+n)*hr(np(26)+n) hr(np(26)+n) = hr(np(111)+n) end do ! n setvar = palloc(111,'TEMP1', 0, 2) endif end