PETSc_dke_ft

Link with the parallel matrix solver package PETSc written in C. 
Adapted to the Saturne computers of the CEA-DRFC

Solve the linear system of equation AX=B

 INPUTS:

   - PETScparam: Structure for controling PETSc program. 
    - A: Matrix [n,n]
   - B: Vector for solving the linear system of equations AX=B [n,1]  
   - display_mode: display some results [1,1]
   - dkepath: paths structure

 OUTPUTS:

   - X: Solution of the linear system of equation [n,1]
   - flag: execution flag (see the MUMPS documentation) [1,1]
   - text: execution documentation (see the HYPRE documentation)

By Joan Decker (CEA-DRFC, joan.decker@cea.fr) ,Yves Peysson (CEA-DRFC,yves.peysson@cea.fr), Frederic Teisseire (CEA-DRFC)

0001 function [X,flag,text] = PETSc_dke_ft(PETScparam,A,B,display_mode,dkepath)
0002 %
0003 %Link with the parallel matrix solver package PETSc written in C.
0004 %Adapted to the Saturne computers of the CEA-DRFC
0005 %
0006 %Solve the linear system of equation AX=B
0007 %
0008 % INPUTS:
0009 %
0010 %   - PETScparam: Structure for controling PETSc program.
0011 %    - A: Matrix [n,n]
0012 %   - B: Vector for solving the linear system of equations AX=B [n,1]
0013 %   - display_mode: display some results [1,1]
0014 %   - dkepath: paths structure
0015 %
0016 % OUTPUTS:
0017 %
0018 %   - X: Solution of the linear system of equation [n,1]
0019 %   - flag: execution flag (see the MUMPS documentation) [1,1]
0020 %   - text: execution documentation (see the HYPRE documentation)
0021 %
0022 %By Joan Decker (CEA-DRFC, joan.decker@cea.fr) ,Yves Peysson (CEA-DRFC,yves.peysson@cea.fr), Frederic Teisseire (CEA-DRFC)
0023 %
0024 if display_mode >= 1,
0025     info_dke_yp(2,'Matrix inversion by PETSc solver - parallel or sequential processing -')
0026 end
0027 %
0028 if nargin < 5,error('Wrong number of input arguments in PETSc_dke_ft !');end
0029 %
0030 flag = 0;text = '';
0031 %
0032 initial_dir = cd;
0033 if ~exist(dkepath.temppath);
0034     mkdir(dkepath.temppath);%Create calculation directory
0035 end
0036 cd(dkepath.temppath);%Move to calculation directory
0037 %
0038 %Write matrix A and vector B in PETSc binary file
0039 %
0040 MAT_FILE_COOKIE = 1211216;
0041 VEC_FILE_COOKIE = 1211214;
0042 %
0043 %Write matrix A
0044 %
0045 fid=fopen('PETSc_in_A.bin','w','ieee-be');
0046 [M,N]=size(A);
0047 fwrite(fid,MAT_FILE_COOKIE,'int');
0048 fwrite(fid,M,'int');
0049 fwrite(fid,N,'int');
0050 fwrite(fid,nnz(A),'int');
0051 [ii,jj,aa]=find(A');
0052 nnzrow=zeros(size(A,1), 1);% nonzeros in each row
0053 nnzrow(:)=full(sum(spones(A'))');
0054 fwrite(fid,nnzrow,'int');
0055 fwrite(fid,ii-1,'int');
0056 fwrite(fid,aa,'float64');
0057 %
0058 %Write vector B
0059 %
0060 fid=fopen('PETSc_in_b.bin','w','ieee-be');  
0061 nv=size(B,2);
0062 for it=1:nv
0063    fwrite(fid,VEC_FILE_COOKIE,'int');
0064    fwrite(fid,size(B,1),'int');
0065    fwrite(fid,full(B(:,it)),'float64');
0066 end
0067 fclose(fid);
0068 %
0069 mattype = PETScparam.matrixtype;
0070 kspmethod = PETScparam.kspmethod;
0071 pcmethod = PETScparam.pcmethod;
0072 %
0073 % Computer management
0074 %
0075 nproc = nproc_dke_yp(display_mode,PETScparam,dkepath);
0076 %
0077 % Run PETSc on N processors (if nproc > 1)
0078 %     - np: specifies number of processors to run on. See mpirun options.
0079 %     - machinefile: specifies a list of available machines to run on. See mpirun options.
0080 %     - f: location of PETSc binary file of the A matrix
0081 %     - rhs: location of PETSc binary file of the b vector
0082 %     - ksp_type: sets KSP method. see PETSc manual
0083 %     - scratchdir: indicates the /scratch directory to the C program.
0084 %
0085 if nproc > 1,%Parallel processing is invoked
0086     if isfield(dkepath.PETSc,'par'),
0087         if display_mode >= 1,info_dke_yp(2,'Matrix inversion by PETSc solver - parallel processing -');end
0088         [flag,text] = unix([dkepath.PETSc.mpipath,'mpirun -np ',int2str(nproc),' -machinefile .machines ',dkepath.PETSc.par,'/PETSc_DKE_ft -f PETSc_in_A.bin -rhs PETSc_in_b.bin  -scratchdir ',dkepath.temppath,' ',mattype,' ',kspmethod, ' ', pcmethod]);
0089     elseif isfield(dkepath.PETSc,'seq'),
0090         if display_mode >= 1,info_dke_yp(2,'Matrix inversion by PETSc solver - sequential processing enforced, no parallel processing package on this machine ! -');end
0091         [flag,text] = unix([dkepath.PETSc.seq,'/PETSc_DKE_ft -f PETSc_in_A.bin -rhs PETSc_in_b.bin  -scratchdir ',dkepath.temppath,' ',mattype,' ',kspmethod, ' ', pcmethod]);
0092     else
0093         error('PETSc package not correctly named for LUKE: <pakage name+version>_P_<compiler name> or  <pakage name+version>_S_<compiler name>');
0094     end
0095 elseif nproc == 1,%Sequential processing is invoked
0096     if isfield(dkepath.PETSc,'seq'),
0097         if display_mode >= 1,info_dke_yp(2,'Matrix inversion by PETSc solver - sequential processing -');end
0098         [flag,text] = unix([dkepath.PETSc.seq,'/PETSc_DKE_ft -f PETSc_in_A.bin -rhs PETSc_in_b.bin  -scratchdir ',dkepath.temppath,' ',mattype,' ',kspmethod, ' ', pcmethod]);
0099     elseif isfield(dkepath.PETSc,'par'),
0100         if display_mode >= 1,info_dke_yp(2,'Matrix inversion by PETSc solver - parallel processing enforced, no sequential processing package on this machine !  -');end
0101         [flag,text] = unix([dkepath.PETSc.mpipath,'mpirun -np ',int2str(nproc),' -machinefile .machines ',dkepath.PETSc.par,'/PETSc_DKE_ft -f PETSc_in_A.bin -rhs PETSc_in_b.bin  -scratchdir ',dkepath.temppath,' ',mattype,' ',kspmethod, ' ', pcmethod]);
0102     else
0103         error('PETSc package not correctly named for LUKE: <pakage name+version>_P_<compiler name> or  <pakage name+version>_S_<compiler name>');
0104     end
0105 end
0106 %
0107 % Read solution from PETSc binary file
0108 %
0109 if ~flag,
0110    fid = fopen('PETSc_out_X.bin','r','ieee-be');
0111    fread(fid,1,'int');
0112    n = fread(fid,1,'int');
0113    X = fread(fid,n,'float64');
0114    fclose(fid);
0115    X = X(:);
0116 else
0117    error(['Error in PETSc calculation. Flag =',int2str(flag),', Log = ',text]);
0118 end
0119 %
0120 delete('*.bin');%delete all bin files
0121 delete('*.bin.info');%delete all bin files
0122 cd(initial_dir);%return to the initial working directory
0123 
0124   
0125   
0126