program gift implicit none save c----------------------------------------------------------------------- c g i f t is a program generator, which outputs fortran code for c g_auss i_nversion of similar linear systems of equations c by a f_ormal t_reatment considering only the non-zero c elements of the corresponding matrices. c c authors: c original version: h. falk, r.g. wolff 1983 c modified by e. mueller april 1996 c f. x. timmes august 2000 c c notes: mueller added the code for a vectorized version. c fxt removed this vectorization for better performance on c parallel machines, and minimized the number of divides. c for n equations, there are now only n divides. c if you want the vectorized version, contact fxt. c c to use this program you will need to modify the code in c routine setmat to reflect the sparsity pattern of your matrix, c and set the logical size of your matrix a few lines below. c you shouldn't have to do much else. c c the output of running this program is set to gift_test.f. c to use a gift generated routine, it is important to note c that the size of the input matrix is a(n,n+1). the extra c column is where the right hand side is stored on input, c and contains the solution to a*x=b on output. c c---------------------------------------------------------------------- c..declare c..physical dimension of the arrays and matrices integer abignet parameter (abignet = 100) character*1 aux(abignet) character*20 string integer n,np,npri,i,j,ilabel,itlsub,nosub,nosub1,numst, 1 lsize,icount,nblk,jj,ii,k,ione,l,istmt double precision a(abignet,abignet), b(abignet) character*9 ctop character*10 csub data ctop/'gift_test'/ data csub/'gsub_test_'/ c..set the logical size of the matrix n = 10 np = n + 1 c..set pattern of matrix a and of vector b call setmat(a, b, n, n, abignet, abignet) c..write out a map of the matrix npri = min0 (n, 132) write(6,6) n 6 format(2x,'map of matrix n=',i4,/) do i = 1, npri do j = 1, npri aux (j) = '.' if (a(i,j) .ne. 0.) aux(j) = 'x' enddo write(6,8) (aux(j), j=1,npri) 8 format(1x,132a1) enddo c..now start writing the fortran code ilabel = 1 ! do-loop label counter itlsub = 1000000 ! maximum number of do-loops per subroutine nosub1 = 100 ! number of first subroutine numst = 0 ! statement number counter within a do-loop lsize = 1000000 ! maximum number of statements per do-loop icount = 0 ! block number counter within a do-loop nblk = 1000000 ! maximum number of blocks per do-loop nosub = nosub1 ! initialize subroutine number counter string = ctop//'.f' open(unit=1,file=string,status='unknown') write (1,3001) csub, nosub,n ! write subroutine header c..forward elimination do 200 jj = 2, n j = n + 2 - jj write (1,6003) j, j, j do 300 ii = 1, j-1 i = j-1 + 1 - ii if (a(i,j) .eq. 0.) goto 11 if (b(j) .eq. 0.) goto 12 if (b(i) .eq. 0.) b(i) = b(j) if ( icount .ge. nblk .and. numst .gt. lsize) then ilabel = ilabel + 1 ! advance do-loop label if ( mod (ilabel, itlsub) .eq. 0) then write(1,2002) ! write return / end nosub = nosub + 1 ! advance subroutine label write(1,3001) csub, nosub,n ! write subroutine header end if numst = 0 icount = 0 end if icount = icount + 1 write (1,1003) i, j, j numst = numst + 1 if (b(i) .eq. 0.) goto 12 write (1,1004) i, np, i, np, j, np numst = numst + 1 12 continue do 400 k = 1, j-1 if (a(j,k) .eq. 0.) goto 400 if (a(i,k) .eq. 0.) a(i,k) = a(j,k) if (a(i,k) .eq. 0.) goto 400 write (1,1005) i, k, i, k, j, k numst = numst + 1 400 continue 11 continue do k = j, n a(i,k) = 0. enddo 300 continue 200 continue ilabel = ilabel + 1 ! advance do-loop counter if ( mod (ilabel, itlsub) .eq. 0) then write(1,2002) ! write return / end nosub = nosub + 1 ! advance subroutine label write(1,3001) csub, nosub,n ! write subroutine header end if icount = 0 numst = 0 c..backsubstitution ione = 1 write (1,1007) ione,ione,ione,ione, np, ione, np, ione do 600 l = 2, n if ( icount .ge. nblk .and. numst .gt. lsize ) then ilabel = ilabel + 1 ! advance do-loop counter if ( mod (ilabel, itlsub) .eq. 0) then write(1,2002) ! write return / end nosub = nosub + 1 ! advance subroutine label write(1,3001) csub, nosub,n ! write subroutine header end if icount = 0 numst = 0 end if icount = icount + 1 write (1,1008) l, np, l, np numst = numst + 1 istmt = 0 do 700 k = 1, l-1 if (a(l,k) .eq. 0.) goto 700 istmt = istmt + 1 if ( mod(istmt,18) .eq. 0) then ! limits continuation lines write (1,4001) ! writes closing ")" numst = numst + 1 icount = icount + 1 write (1,1008) l, np, l, np ! new statement numst = numst + 1 end if write (1,1009) l, k, k, np ! continuation card numst = numst + 1 700 continue write (1,2001) l numst = numst + 1 600 continue write (1,2002) ! write return / end write (1,1001) ctop ! write master subroutine header do k = nosub1, nosub write (1,3002) csub, k ! write call to sub-subroutine enddo write (1,2002) ! write return / end close(unit=1) print *, ' ' print *, ' total number of do-loops: ', ilabel print *, ' total number of subroutines: ', nosub - nosub1 + 1 print *, ' subroutines with prefix: ',csub print *, ' and subroutine gift have been written to ',string print *, ' last subroutine is ',csub,nosub stop c..format statements 1001 format(' subroutine ',a,'(ab,n1,n2)'/ & ' implicit none'/ & ' integer n1,n2'/ & ' double precision ab(n1,n2)' ) 3001 format(' subroutine ',a, i3, '(ab,n1,n2)'/ & ' implicit none'/ & ' integer n1,n2'/ & ' double precision ab(n1,n2),c,tmp(',i3,')') 3002 format(' call ',a, i3, '(ab,n1,n2)') 1002 format('c'/ ' do ', i4, ' k = 1, np') 6003 format(1x,/,8x,'tmp(',i3,') = 1.0d0/ab(',i3,',',i3,')') 1003 format(8x,'c', 11x, '= ab(',i3,',',i3,') * tmp(',i4,')') 1004 format(8x,'ab(', i3, ',', i3,') = ab(', i3, ',', i3 ,')', & ' - ab(', i3, ',', i3, ') * c') 1005 format(8x,'ab(', i3, ',', i3, ') = ab(', i3, ',', i3, ')', & ' - ab(', i3, ',', i3, ') * c') 1006 format(1x, i4, ' continue') 1007 format(1x,/,1x,/,'c backsubstitution',/, & 8x,'tmp(',i3,') = 1.0d0/ab(',i3,',',i3,')',/, & 8x,'ab(', i3, ',', i3, ') = ab(', i3, ',', i3, ')', & ' * tmp(', i4, ')' ) 1008 format(8x,'ab(', i3, ',', i3, ') = ( ab(', i3, ',', i3, ')' ) 1009 format(5x, '&', 17x, ' - ab(', i3, ',', i3, ')', & ' * ab(', i3, ',', i3, ')' ) 2001 format(5x, '&', 46x, ') *tmp(', i4, ')' ) 4001 format(5x, '&', 46x, ')' ) 2002 format(6x, 'return'/6x, 'end') end subroutine setmat(a, b, nl, n, nlphys, nphys) implicit none save c..this routine sets the nonzero pattern of the dense matrix a, c..whose physical dimension is a(nlphys,nphys) and logical dimension a(nl,n) c..the nonzero pattern of the dense right hand side b, c..whose physical dimension is b(nlphys) and logical dimension b(n) c..declare integer np parameter (np=20000) integer nl,n,nlphys,nphys,i,j,inew,jnew, 1 iloc(np),jloc(np),nzo double precision a(nlphys,nphys),b(nphys),xx,yy c..initialize the matrix and right hand side do j=1,n do i=1,nl a(i,j) = 0.0d0 enddo enddo do i=1,n b(i) = 1.0d0 enddo c..here i set the pattern for a simple tridiagonal c..in a real application you must modify this portion of the program a(1,1) = 1.0d0 a(1,2) = 1.0d0 do j=2,n-1 do i=j-1,j+1 a(j,i) = 1.0d0 enddo enddo a(n,n-1) = 1.0d0 a(n,n) = 1.0d0 return end