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 PetscErrorCode ierr; 44 int m,n = 0,nz,dummy; /* these are fscaned so kept as int */ 45 PetscInt i,col,row,shift = 1,sizes[3],nsizes; 46 PetscScalar val; 47 PetscReal res_norm; 48 FILE *Afile,*bfile,*ufile; 49 PetscViewer view; 50 PetscBool flg_A,flg_b,flg_u,flg; 51 PetscMPIInt size; 52 53 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 54 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 55 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 56 57 /* Read in matrix, rhs and exact solution from ascii files */ 58 ierr = PetscOptionsGetString(NULL,NULL,"-Ain",Ain,sizeof(Ain),&flg_A);CHKERRQ(ierr); 59 ierr = PetscOptionsHasName(NULL,NULL,"-noshift",&flg);CHKERRQ(ierr); 60 if (flg) shift = 0; 61 if (flg_A) { 62 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");CHKERRQ(ierr); 63 ierr = PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);CHKERRQ(ierr); 64 nsizes = 3; 65 ierr = PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg);CHKERRQ(ierr); 66 if (flg) { 67 if (nsizes != 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must pass in three m,n,nz as arguments for -nosizesinfile"); 68 m = sizes[0]; 69 n = sizes[1]; 70 nz = sizes[2]; 71 } else { 72 if (fscanf(Afile,"%d %d %d\n",&m,&n,&nz) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 73 } 74 ierr = PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);CHKERRQ(ierr); 75 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example\n"); 76 ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr); 77 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 78 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 79 ierr = MatSeqAIJSetPreallocation(A,nz/m,NULL);CHKERRQ(ierr); 80 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 81 82 for (i=0; i<nz; i++) { 83 if (fscanf(Afile,"%d %d %le\n",&row,&col,(double*)&val) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 84 row -= shift; col -= shift; /* set index set starts at 0 */ 85 ierr = MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);CHKERRQ(ierr); 86 } 87 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 88 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 89 fclose(Afile); 90 } 91 92 ierr = PetscOptionsGetString(NULL,NULL,"-rhs",rhs,sizeof(rhs),&flg_b);CHKERRQ(ierr); 93 if (flg_b) { 94 ierr = VecCreate(PETSC_COMM_SELF,&b);CHKERRQ(ierr); 95 ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr); 96 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 97 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");CHKERRQ(ierr); 98 ierr = PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);CHKERRQ(ierr); 99 for (i=0; i<n; i++) { 100 if (fscanf(bfile,"%d %le\n",&dummy,(double*)&val) != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 101 ierr = VecSetValues(b,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr); 102 } 103 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 104 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); 105 fclose(bfile); 106 } 107 108 ierr = PetscOptionsGetString(NULL,NULL,"-solu",solu,sizeof(solu),&flg_u);CHKERRQ(ierr); 109 if (flg_u) { 110 ierr = VecCreate(PETSC_COMM_SELF,&u);CHKERRQ(ierr); 111 ierr = VecSetSizes(u,PETSC_DECIDE,n);CHKERRQ(ierr); 112 ierr = VecSetFromOptions(u);CHKERRQ(ierr); 113 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");CHKERRQ(ierr); 114 ierr = PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);CHKERRQ(ierr); 115 for (i=0; i<n; i++) { 116 if (fscanf(ufile,"%d %le\n",&dummy,(double*)&val) != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 117 ierr = VecSetValues(u,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr); 118 } 119 ierr = VecAssemblyBegin(u);CHKERRQ(ierr); 120 ierr = VecAssemblyEnd(u);CHKERRQ(ierr); 121 fclose(ufile); 122 } 123 124 /* Write matrix, rhs and exact solution in Petsc binary file */ 125 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n");CHKERRQ(ierr); 126 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);CHKERRQ(ierr); 127 ierr = MatView(A,view);CHKERRQ(ierr); 128 129 if (flg_b) { /* Write rhs in Petsc binary file */ 130 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n");CHKERRQ(ierr); 131 ierr = VecView(b,view);CHKERRQ(ierr); 132 } 133 if (flg_u) { 134 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n");CHKERRQ(ierr); 135 ierr = VecView(u,view);CHKERRQ(ierr); 136 } 137 138 /* Check accuracy of the data */ 139 if (flg_A & flg_b & flg_u) { 140 ierr = VecDuplicate(u,&u_tmp);CHKERRQ(ierr); 141 ierr = MatMult(A,u,u_tmp);CHKERRQ(ierr); 142 ierr = VecAXPY(u_tmp,-1.0,b);CHKERRQ(ierr); 143 ierr = VecNorm(u_tmp,NORM_2,&res_norm);CHKERRQ(ierr); 144 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);CHKERRQ(ierr); 145 ierr = VecDestroy(&u_tmp);CHKERRQ(ierr); 146 } 147 148 ierr = MatDestroy(&A);CHKERRQ(ierr); 149 if (flg_b) {ierr = VecDestroy(&b);CHKERRQ(ierr);} 150 if (flg_u) {ierr = VecDestroy(&u);CHKERRQ(ierr);} 151 ierr = PetscViewerDestroy(&view);CHKERRQ(ierr); 152 ierr = PetscFinalize(); 153 return ierr; 154 } 155 156 /*TEST 157 158 build: 159 requires: !defined(PETSC_USE_64BIT_INDICES) double !complex 160 161 test: 162 requires: datafilespath 163 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat 164 165 TEST*/ 166