1 2 static char help[] = "Reads in a matrix in ASCII MATLAB format (I,J,A), read in vectors rhs and exact_solu in ASCII format.\n\ 3 Writes them using the PETSc sparse format.\n\ 4 Note: I and J start at 1, not 0, use -noshift if indices in file start with zero!\n\ 5 Input parameters are:\n\ 6 -Ain <filename> : input matrix in ascii format\n\ 7 -rhs <filename> : input rhs in ascii format\n\ 8 -solu <filename> : input true solution in ascii format\n\\n"; 9 10 /* 11 Example: ./ex78 -Ain Ain -rhs rhs -solu solu -noshift -mat_view 12 with the datafiles in the following format: 13 Ain (I and J start at 0): 14 ------------------------ 15 3 3 6 16 0 0 1.0 17 0 1 2.0 18 1 0 3.0 19 1 1 4.0 20 1 2 5.0 21 2 2 6.0 22 23 rhs 24 --- 25 0 3.0 26 1 12.0 27 2 6.0 28 29 solu 30 ---- 31 0 1.0 32 0 1.0 33 0 1.0 34 */ 35 36 #include <petscmat.h> 37 38 int main(int argc,char **args) 39 { 40 Mat A = NULL; 41 Vec b,u = NULL,u_tmp; 42 char Ain[PETSC_MAX_PATH_LEN],rhs[PETSC_MAX_PATH_LEN],solu[PETSC_MAX_PATH_LEN]; 43 int m,n = 0,nz,dummy; /* these are fscaned so kept as int */ 44 PetscInt i,col,row,shift = 1,sizes[3],nsizes; 45 PetscScalar val; 46 PetscReal res_norm; 47 FILE *Afile,*bfile,*ufile; 48 PetscViewer view; 49 PetscBool flg_A,flg_b,flg_u,flg; 50 PetscMPIInt size; 51 52 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 53 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 54 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 55 56 /* Read in matrix, rhs and exact solution from ascii files */ 57 PetscCall(PetscOptionsGetString(NULL,NULL,"-Ain",Ain,sizeof(Ain),&flg_A)); 58 PetscCall(PetscOptionsHasName(NULL,NULL,"-noshift",&flg)); 59 if (flg) shift = 0; 60 if (flg_A) { 61 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n")); 62 PetscCall(PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile)); 63 nsizes = 3; 64 PetscCall(PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg)); 65 if (flg) { 66 PetscCheckFalse(nsizes != 3,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must pass in three m,n,nz as arguments for -nosizesinfile"); 67 m = sizes[0]; 68 n = sizes[1]; 69 nz = sizes[2]; 70 } else { 71 PetscCheckFalse(fscanf(Afile,"%d %d %d\n",&m,&n,&nz) != 3,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 72 } 73 PetscCall(PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz)); 74 PetscCheckFalse(m != n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example"); 75 PetscCall(MatCreate(PETSC_COMM_SELF,&A)); 76 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n)); 77 PetscCall(MatSetFromOptions(A)); 78 PetscCall(MatSeqAIJSetPreallocation(A,nz/m,NULL)); 79 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 80 81 for (i=0; i<nz; i++) { 82 PetscCheckFalse(fscanf(Afile,"%d %d %le\n",&row,&col,(double*)&val) != 3,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 83 row -= shift; col -= shift; /* set index set starts at 0 */ 84 PetscCall(MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES)); 85 } 86 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 87 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 88 fclose(Afile); 89 } 90 91 PetscCall(PetscOptionsGetString(NULL,NULL,"-rhs",rhs,sizeof(rhs),&flg_b)); 92 if (flg_b) { 93 PetscCall(VecCreate(PETSC_COMM_SELF,&b)); 94 PetscCall(VecSetSizes(b,PETSC_DECIDE,n)); 95 PetscCall(VecSetFromOptions(b)); 96 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n")); 97 PetscCall(PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile)); 98 for (i=0; i<n; i++) { 99 PetscCheckFalse(fscanf(bfile,"%d %le\n",&dummy,(double*)&val) != 2,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 100 PetscCall(VecSetValues(b,1,&i,&val,INSERT_VALUES)); 101 } 102 PetscCall(VecAssemblyBegin(b)); 103 PetscCall(VecAssemblyEnd(b)); 104 fclose(bfile); 105 } 106 107 PetscCall(PetscOptionsGetString(NULL,NULL,"-solu",solu,sizeof(solu),&flg_u)); 108 if (flg_u) { 109 PetscCall(VecCreate(PETSC_COMM_SELF,&u)); 110 PetscCall(VecSetSizes(u,PETSC_DECIDE,n)); 111 PetscCall(VecSetFromOptions(u)); 112 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n")); 113 PetscCall(PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile)); 114 for (i=0; i<n; i++) { 115 PetscCheckFalse(fscanf(ufile,"%d %le\n",&dummy,(double*)&val) != 2,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 116 PetscCall(VecSetValues(u,1,&i,&val,INSERT_VALUES)); 117 } 118 PetscCall(VecAssemblyBegin(u)); 119 PetscCall(VecAssemblyEnd(u)); 120 fclose(ufile); 121 } 122 123 /* Write matrix, rhs and exact solution in Petsc binary file */ 124 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n")); 125 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view)); 126 PetscCall(MatView(A,view)); 127 128 if (flg_b) { /* Write rhs in Petsc binary file */ 129 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n")); 130 PetscCall(VecView(b,view)); 131 } 132 if (flg_u) { 133 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n")); 134 PetscCall(VecView(u,view)); 135 } 136 137 /* Check accuracy of the data */ 138 if (flg_A & flg_b & flg_u) { 139 PetscCall(VecDuplicate(u,&u_tmp)); 140 PetscCall(MatMult(A,u,u_tmp)); 141 PetscCall(VecAXPY(u_tmp,-1.0,b)); 142 PetscCall(VecNorm(u_tmp,NORM_2,&res_norm)); 143 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm)); 144 PetscCall(VecDestroy(&u_tmp)); 145 } 146 147 PetscCall(MatDestroy(&A)); 148 if (flg_b) PetscCall(VecDestroy(&b)); 149 if (flg_u) PetscCall(VecDestroy(&u)); 150 PetscCall(PetscViewerDestroy(&view)); 151 PetscCall(PetscFinalize()); 152 return 0; 153 } 154 155 /*TEST 156 157 build: 158 requires: !defined(PETSC_USE_64BIT_INDICES) double !complex 159 160 test: 161 requires: datafilespath 162 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat 163 164 TEST*/ 165