home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / skyline / example / dskstest.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  7.7 KB  |  319 lines
  1.       program main
  2. c
  3. c     This is an example program to test the symmetric skyline solver.
  4. c
  5. c     Input:
  6. c            n        dimension of the problem
  7. c            nb        average band width desired when generating data
  8. c            nv         the std. dev. wanted for the band width
  9. c            igen       = 0      the matrix is generated randomly
  10. c                       > or = 1 the position of the diagonal elements
  11. c                                must be given
  12.  
  13.       implicit double precision (a-h,o-z)
  14.  
  15.       integer          id(10001)
  16.       double precision a(1000000),ao(1000000),b(10000)
  17.  
  18. 100   print *, ' '
  19.       print *, 'Enter n, average bandwidth, bandwidth variance:'     
  20.       read( 5, *, end=999 ) n, nb, nv
  21.       print *, 'Storage Scheme: Enter  1  for NORMAL Jennings '
  22.       print *, '                       2  for REVERSE Jennings:'
  23.       read( 5, *, end=999 ) ityp
  24.       print *, 'Matrix to be generated randomly: '
  25.       print *, '       Enter  0  for whole matrix generated randomly'
  26.       print *, ' '
  27.       print *, '              1  for diagonal elements given,' 
  28.       print *, '                     matrix generated randomly'
  29.       print *, ' '
  30.       print *, '              2  for diagonal elements given,'
  31.       print *, '                     matrix element given :'
  32.       read( 5, *, end=999 ) igen
  33.  
  34. c     Use as many CPU as available
  35.       ncpu = mp_numthreads() 
  36.  
  37. c     Use Jennings Method because reverse Jennings is NOT
  38. c     available yet.
  39.       print *, ' '
  40.       if( ityp .EQ. 1 ) then
  41.       print *,'main:  jenning\'s method '
  42.       else
  43.       print *,'main:  reverse jenning\'s method '
  44.       end if
  45.  
  46.       print *,'       num cpu in system = num cpu used: ', ncpu
  47.  
  48.       call driver( n, nb, nv, igen, ityp, a, ao, b, id )
  49.  
  50.       go to 100
  51.  
  52. 999   stop
  53.       end
  54.  
  55.       subroutine driver( n, nb, nv, igen, ityp, a, ao, b, id )
  56.  
  57.       implicit double precision (a-h,o-z)
  58.  
  59.       integer           n, nb, nv, igen, ityp
  60.       integer           id(*)
  61.       double precision  a(*),ao(*),b(*)
  62.  
  63.       integer           am
  64.       integer           iops
  65.       real              etime,t1,et,wt,tarray(2),et1,second
  66.       real*8            rand
  67.       double precision  tol
  68.       double precision  zero, one, two, half, onem
  69.       double precision  tmp
  70.  
  71.       parameter         ( zero = 0.0d0 )
  72.       parameter         ( one  = 1.0d0 )
  73.       parameter         ( two  = 2.0d0 )
  74.       parameter         ( half = 5.0d-1 )
  75.       parameter         ( oneh = 1.0d2 )
  76.       parameter         ( onem = 1.0d6 )
  77.  
  78.       integer           idmach
  79.       double precision  dmach
  80.  
  81.       intsize = idmach( 'I' )
  82.       iflsize = idmach( 'F' )
  83.       tol     = dmach( 'E' )
  84.  
  85.       if(igen .eq. 0) then
  86.         iseed=10001
  87.         call srand(iseed)
  88.         do i=1,n
  89.            a(i) = rand()
  90.         enddo
  91.  
  92.         id(1)=1
  93. c    if( ityp .EQ. 1 ) then
  94.            do i=2,n
  95.               col= abs(two*(a(i)-half)*float(nv)+float(nb)+half)
  96.               ilen=min0(int(col),i)
  97.               ilen=max0(ilen,1)
  98.               id(i)=id(i-1)+ilen
  99.            enddo
  100. c    else
  101. c       id(2) = 2
  102. c           do i=2,n
  103. c              col= abs(two*(a(i)-half)*float(nv)+float(nb)+half)
  104. c              ilen=min0(int(col),i)
  105. c              ilen=max0(ilen,1)
  106. c              id(i+1)=id(i)+ilen+1
  107. c           enddo
  108. c        end if
  109.  
  110.       elseif(igen .ge. 1) then
  111.            read(5,*) (id(i),i=1,n)
  112.       endif
  113.  
  114.       am=id(1)
  115.       hm=one/float(id(1))
  116.       mx=id(1)
  117.  
  118.       do i=2,n
  119.         j=id(i)-id(i-1)
  120.         mx=max(j,mx)
  121.         am=am+j
  122.         hm=hm+one/float(j)
  123.       enddo
  124.  
  125.       print *, ' '
  126.       print*,'       ndof: ',n
  127.       ip=100.*float(mx)/float(n)+.5
  128.       print*,'       max band: ',mx,' or ',ip,' % '
  129.  
  130.       am=am/float(n)
  131.       iam=am
  132.       ip=100.*am/float(n)+.5
  133.       print*,'       arithmetic mean band: ',iam,' or ',ip,' % '
  134.  
  135.       hm=hm/float(n)
  136.       hm=one/hm
  137.       ihm=hm
  138.       ip=100.*hm/float(n)+.5
  139.       print*,'         harmonic mean band: ',ihm,' or ',ip,' % '
  140.  
  141.       lena=id(n)
  142.       i=float(8*lena)/onem+.5
  143.       print*,'       num of terms: ',lena,' megabytes: ',i
  144.  
  145.       if(igen .le. 1) then
  146.         do i=1,lena
  147.            a(i) = rand()
  148.         enddo
  149.         do i=1,n
  150.           a(id(i))=(a(id(i))+half)*oneh
  151.         enddo
  152.       elseif(igen .eq. 2) then
  153.         read(5,*) (a(i),i=1,lena)
  154.       endif
  155.  
  156.       do i=1,lena
  157.         ao(i)=a(i)
  158.       enddo
  159.       do i=1,n
  160.         b(i)=one
  161.       enddo
  162.  
  163.       if( ityp .EQ. 2 ) then
  164.       id( n+1 ) = id( n ) + 1
  165.       do j = n, 2, -1
  166.          ibgn = id( j-1 ) + 1
  167.          iend = id( j )
  168.          id( j ) = ibgn
  169.          imax = ( iend - ibgn + 1 )/2
  170.          do i = 1, imax
  171.         tmp = a( iend )
  172.         a( iend ) = a( ibgn )
  173.         a( ibgn ) = tmp
  174.         ibgn = ibgn + 1
  175.         iend = iend - 1
  176.          end do
  177.       end do
  178.       end if
  179.  
  180.          et1=second()
  181.          t1=etime(tarray)
  182.          call dskydc(a,n,id,info,ityp)
  183.          if(info .gt. 0)  then
  184.         print *, '**** INFO = ', info
  185.         return
  186.      end if
  187.      if(info .LT. 0) return
  188.  
  189.          et=etime(tarray) - t1
  190.          wt=second()-et1
  191.  
  192.  
  193.  
  194.       call djesol(a,n,id,b)
  195.       call error(ao,n,id,b)
  196.  
  197.       big=zero
  198.       do i=1,n
  199.         err=abs(b(i)-one)
  200.         if(err .gt. big) then
  201.           big=err
  202.         endif
  203.       enddo
  204.       print*,'       error= ',big
  205.       
  206.       if(nv .eq. 0) then
  207.         xb=nb
  208.         xn=n
  209.         ops=xb*(xb-1.)*(3.*xn-2.*xb+1.)/3.
  210.         ops=ops/onem+.5
  211.       else
  212.         iops=0
  213.         call counter(id,n,iops)
  214.         ops=float(iops)/onem+.5
  215.       endif
  216.       if(et .gt. zero) flops=ops/et
  217.       if(wt .gt. zero) then
  218.          wflops=ops/wt
  219.          ratio=wt/et
  220.       endif
  221.  
  222.       write(6,48) ops
  223. 48    format(/,' million fp ops = ',f8.2)
  224.       write(6,49) et,wt
  225. 49    format(' cpu time  = ',f6.2,' wall time  = ',f6.2)
  226.       write(6,50) flops,wflops,ratio
  227. 50    format(' cpu mflops= ',f6.2,' wall mflops= ',f6.2,' ratio= ',f6.2)
  228.  
  229.       return
  230.       end
  231.  
  232.       SUBROUTINE COUNTER(LC,NEQ,ICOUNT)
  233. C---------------------------------------------- 
  234. C      ---ACTIVE COLUMN EQUATION SOLVER  AX=B-- 
  235. C          KK=1 LDL FACTORAZATION ONLY
  236. C          KK=2 FORWARD REDUCTIION ONLY 
  237. C          KK=3 BACKSUBSTITION ONLY 
  238. C---------------------------------------------- 
  239.       DIMENSION LC(NEQ)
  240.       DOUBLE PRECISION S,T,FF
  241.       ICOUNT = 0
  242.  50   IF(NEQ.EQ.1) RETURN 
  243. C------LDL FACTORAZATION--------------------- 
  244.       DO 500 J=2,NEQ
  245.       JH=LC(J)-LC(J-1)
  246.       IF(JH.EQ.1) GO TO 500 
  247.       K=J-JH+1
  248. C-----FORM U(I,J) - TOP OF COLUMN DOWN TO DIAGONAL
  249.       I=K 
  250.       LCIM1=0
  251.       IF(I-1 .GT. 0) LCIM1=LC(I-1)
  252.  100  NT=MIN0(JH-J+I,LC(I)-LCIM1)-1 
  253.       IF(NT) 300,300,150
  254.  150  NS=LC(I)-NT 
  255.       NE=LC(I)-1
  256.       IJ=LC(J)-J+I
  257.       IC=IJ-LC(I) 
  258. C**      S=0.0D0 
  259.       IF(I.EQ.J) GO TO 400
  260.       ICOUNT = ICOUNT + 2 * NT 
  261. C**      CALL DOTP(A(NS),A(NS+IC),S,NT)
  262. C**      A(IJ)=A(IJ)-S 
  263.  300  I=I+1 
  264.       GO TO 100 
  265. C-----FORM L(I,J) AND U(I,I) ------------------ 
  266.  400  DO 450 N=NS,NE
  267.       ND=LC(K)
  268.       K=K+1 
  269. C**      T=A(N)
  270. C**      A(N)=A(N)/A(ND) 
  271.       ICOUNT = ICOUNT + 3
  272. C** 450  S=S+A(N)*T
  273.  450  CONTINUE
  274. C**      A(IJ)=A(IJ)-S 
  275.       ICOUNT = ICOUNT + 1
  276.  500  CONTINUE
  277.       RETURN
  278.       END 
  279.  
  280.       SUBROUTINE ERROR(A,N,ID,X)
  281.       POINTER (PLHT,LHT), (PZ,Z)
  282.       INTEGER  N,ID(1),LHT(1),M
  283.       DOUBLE PRECISION A(1),X(1),Z(1),S,TOL
  284.       INTEGER          IDMACH
  285.       DOUBLE PRECISION DMACH
  286.  
  287.       INTSIZE = IDMACH( 'I' )
  288.       IFLSIZE = IDMACH( 'F' )
  289.       TOL     = DMACH( 'E' )
  290.  
  291.       PLHT=MALLOC(N*INTSIZE)
  292.       PZ=MALLOC(N*IFLSIZE)
  293.       LHT(1)=0
  294.       DO 5 I=2,N
  295.         LHT(I)=ID(I)-ID(I-1)-1
  296.     5 CONTINUE
  297.       DO 10 K=1,N
  298.         Z(K)=X(K)*A(ID(K))
  299.    10 CONTINUE
  300.       DO 40 K=2,N
  301.         M=ID(K-1)+1
  302.         II=0
  303.         S=Z(K)
  304.         DO 30 I=K-LHT(K),K-1
  305.           Z(I)=Z(I)+X(K)*A(M+II)
  306.           S=S+X(I)*A(M+II)
  307.           II=II+1
  308.    30   CONTINUE
  309.         Z(K)=S
  310.    40 CONTINUE
  311.       DO 50 I=1,N
  312.         X(I) = Z(I)
  313.    50 CONTINUE
  314.       CALL FREE(PLHT)
  315.       CALL FREE(PZ)
  316.       RETURN
  317.       END
  318.  
  319.