mumpsdke

Fortran source

0001 C This file must be used with DMUMPS VERSION 4.6.1
0002 C
0003 CCCCCCCCCCC ADAPTED FOR THE ELECTRON DRIFT KINETIC SOLVER DKE CCCCCCCCCCCCC 
0004 CC By Y.PEYSSON (yves.peysson@cea.fr) and J. DECKER (joan.decker@cea.fr) CC
0005 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0006 C
0007       PROGRAM DMUMPS_DKE
0008       IMPLICIT NONE
0009       INCLUDE 'mpif.h'
0010       INCLUDE 'dmumps_struc.h'
0011       TYPE (DMUMPS_STRUC) mumps_par
0012       INTEGER IERR, I
0013 C
0014 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FOR DKE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
0015 C
0016       double precision:: dimvalu
0017       integer:: rw,rc,lengthin,lengthout
0018       real*8 ZZ
0019 C
0020 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0021 C
0022       
0023       CALL MPI_INIT(IERR)
0024 C Define a communicator for the package.
0025       mumps_par%COMM = MPI_COMM_WORLD
0026 C  Initialize an instance of the package
0027 C  for L U factorization (sym = 0, with working host)
0028       mumps_par%JOB = -1
0029       mumps_par%SYM = 0
0030       mumps_par%PAR = 1
0031       CALL DMUMPS(mumps_par)
0032 C  Define problem on the host (processor 0)
0033       IF ( mumps_par%MYID .eq. 0 ) THEN
0034 
0035 C
0036 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FOR DKE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
0037 C
0038 
0039       inquire(iolength=lengthin) dimvalu
0040       
0041 C  Read vector size
0042 
0043       open(unit=77,file='MUMPS_N_in.bin',access='direct',
0044      *form='unformatted',status='old',recl=lengthin,iostat=rc,err=300)      
0045       read(unit=77,rec=1,iostat=rw) ZZ
0046       mumps_par%N = ZZ
0047       close(77)
0048 C      write(*,*) ZZ,mumps_par%N
0049       
0050 C  Read number of non-zero elements in matrix A 
0051 
0052       open(unit=77,file='MUMPS_NA_in.bin',access='direct',
0053      *form='unformatted',status='old',recl=lengthin,iostat=rc,err=301)      
0054       read(unit=77,rec=1,iostat=rw) ZZ
0055       mumps_par%NZ = ZZ
0056       close(77)
0057       
0058 C  Memory allocation
0059 
0060       ALLOCATE( mumps_par%IRN ( mumps_par%NZ ) )
0061       ALLOCATE( mumps_par%JCN ( mumps_par%NZ ) )
0062       ALLOCATE( mumps_par%A( mumps_par%NZ ) )
0063       ALLOCATE( mumps_par%RHS ( mumps_par%N  ) )      
0064 
0065 C  Read non-zero elements of matrix A      
0066       
0067       open(unit=77,file='MUMPS_AI_in.bin',access='direct',
0068      *form='unformatted',status='old',recl=lengthin,iostat=rc,err=302) 
0069       open(unit=78,file='MUMPS_AJ_in.bin',access='direct',
0070      *form='unformatted',status='old',recl=lengthin,iostat=rc,err=303) 
0071       open(unit=79,file='MUMPS_AV_in.bin',access='direct',
0072      *form='unformatted',status='old',recl=lengthin,iostat=rc,err=304) 
0073       
0074       DO I = 1, mumps_par%NZ
0075          read(unit=77,rec=I,iostat=rw) ZZ
0076          mumps_par%IRN(I) = ZZ
0077          read(unit=78,rec=I,iostat=rw) ZZ
0078          mumps_par%JCN(I) = ZZ
0079          read(unit=79,rec=I,iostat=rw) ZZ
0080          mumps_par%A(I) = ZZ
0081       END DO
0082       
0083       close(77)
0084       close(78)
0085       close(79)
0086 
0087 C  Read vector B
0088 
0089       open(unit=77,file='MUMPS_B_in.bin',access='direct',
0090      *form='unformatted',status='old',recl=lengthin,iostat=rc,err=305) 
0091 
0092       DO I = 1, mumps_par%N
0093          read(unit=77,rec=I,iostat=rw) ZZ
0094          mumps_par%RHS(I) = ZZ  
0095       END DO
0096 
0097       close(77)
0098 
0099       END IF
0100 C  Call package for solution
0101       mumps_par%JOB = 6
0102       CALL DMUMPS(mumps_par)
0103 C  Solution has been assembled on the host
0104       IF ( mumps_par%MYID .eq. 0 ) THEN
0105 C
0106 CCCCCCCCCCCCCCCCCCCCCCCCCCCC FOR DKE OUTPUT CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
0107 C
0108       inquire(iolength=lengthout) dimvalu
0109 
0110 C  Write vector X
0111 
0112       open(unit=77,file='MUMPS_X_out.bin',access='direct',
0113      *form='unformatted',status='replace',
0114      *recl=lengthout,iostat=rc,err=400)
0115       
0116       
0117       do i=1,mumps_par%N
0118         write(unit=77,rec=I,iostat=rw) mumps_par%RHS(I)
0119       enddo
0120 
0121       close(77)      
0122       
0123       GOTO 1000
0124       
0125 C
0126 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0127 C
0128       END IF
0129 C  Deallocate user data
0130       
0131 C
0132 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FOR DKE   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0133 C
0134 
0135       
0136   300 continue
0137       print*,'PROBLEM IN READING INPUT FILE MUMPS_N_in.bin'
0138       GOTO 1000
0139   301 continue
0140       print*,'PROBLEM IN READING INPUT FILE MUMPS_NA_in.bin'
0141       GOTO 1000
0142   302 continue
0143       print*,'PROBLEM IN READING INPUT FILE MUMPS_AI_in.bin'
0144       GOTO 1000
0145   303 continue
0146       print*,'PROBLEM IN READING INPUT FILE MUMPS_AJ_in.bin'
0147       GOTO 1000
0148   304 continue
0149       print*,'PROBLEM IN READING INPUT FILE MUMPS_AV_in.bin'
0150       GOTO 1000
0151   305 continue
0152       print*,'PROBLEM IN READING INPUT FILE MUMPS_B_in.bin'
0153       GOTO 1000
0154       
0155   400 continue
0156       print*,
0157      *'PROBLEM IN WRITING OUTPUT FILE MUMPS_X_out.bin'
0158       GOTO 1000
0159 
0160 C
0161 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0162 C
0163 
0164  1000 IF ( mumps_par%MYID .eq. 0 )THEN
0165         DEALLOCATE( mumps_par%IRN )
0166         DEALLOCATE( mumps_par%JCN )
0167         DEALLOCATE( mumps_par%A   )
0168         DEALLOCATE( mumps_par%RHS )
0169       END IF
0170 C  Destroy the instance (deallocate internal data structures)
0171       mumps_par%JOB = -2
0172       CALL DMUMPS(mumps_par)
0173       CALL MPI_FINALIZE(IERR)
0174 
0175       STOP
0176       END
0177