1 2 static char help[] = "Tests the use of interface functions for MATIS matrices.\n\ 3 This example tests: MatZeroRows(), MatZeroRowsLocal(), MatView(), MatDuplicate(),\n\ 4 MatCopy(), MatCreateSubMatrix(), MatGetLocalSubMatrix(), MatAXPY(), MatShift()\n\ 5 MatDiagonalSet(), MatTranspose() and MatPtAP(). It also tests some\n\ 6 conversion routines.\n"; 7 8 #include <petscmat.h> 9 10 PetscErrorCode TestMatZeroRows(Mat,Mat,PetscBool,IS,PetscScalar); 11 PetscErrorCode CheckMat(Mat,Mat,PetscBool,const char*); 12 13 int main(int argc,char **args) 14 { 15 Mat A,B,A2,B2,T; 16 Mat Aee,Aeo,Aoe,Aoo; 17 Mat *mats; 18 Vec x,y; 19 MatInfo info; 20 ISLocalToGlobalMapping cmap,rmap; 21 IS is,is2,reven,rodd,ceven,codd; 22 IS *rows,*cols; 23 MatType lmtype; 24 PetscScalar diag = 2.; 25 PetscInt n,m,i; 26 PetscInt rst,ren,cst,cen,nr,nc; 27 PetscMPIInt rank,size; 28 PetscBool testT,squaretest,isaij; 29 PetscBool permute = PETSC_FALSE; 30 PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE; 31 PetscErrorCode ierr; 32 33 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 34 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 35 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 36 m = n = 2*size; 37 ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);CHKERRQ(ierr); 38 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 40 if (size > 1 && m < 4) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 4 for parallel runs"); 41 if (size == 1 && m < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 2 for uniprocessor runs"); 42 if (n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of cols should be larger or equal 2"); 43 if (symmetric) m = n = PetscMax(m,n); 44 45 /* create a MATIS matrix */ 46 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 47 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 48 ierr = MatSetType(A,MATIS);CHKERRQ(ierr); 49 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 50 /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines */ 51 ierr = ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is);CHKERRQ(ierr); 52 ierr = ISLocalToGlobalMappingCreateIS(is,&cmap);CHKERRQ(ierr); 53 ierr = ISDestroy(&is);CHKERRQ(ierr); 54 55 ierr = PetscOptionsGetBool(NULL,NULL,"-permmap",&permute,NULL);CHKERRQ(ierr); 56 ierr = PetscOptionsGetBool(NULL,NULL,"-diffmap",&diffmap,NULL);CHKERRQ(ierr); 57 if (!symmetric && (diffmap || m != n)) { 58 59 ierr = ISCreateStride(PETSC_COMM_WORLD,m,permute ? m -1 : 0,permute ? -1 : 1,&is);CHKERRQ(ierr); 60 ierr = ISLocalToGlobalMappingCreateIS(is,&rmap);CHKERRQ(ierr); 61 ierr = ISDestroy(&is);CHKERRQ(ierr); 62 if (m==n && !permute) squaretest = PETSC_TRUE; 63 else squaretest = PETSC_FALSE; 64 } else { 65 ierr = PetscObjectReference((PetscObject)cmap);CHKERRQ(ierr); 66 rmap = cmap; 67 squaretest = PETSC_TRUE; 68 } 69 ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr); 70 ierr = MatISStoreL2L(A,PETSC_FALSE);CHKERRQ(ierr); 71 ierr = MatISSetPreallocation(A,3,NULL,0,NULL);CHKERRQ(ierr); 72 for (i=0; i<m; i++) { 73 PetscScalar v[3]; 74 PetscInt cols[3]; 75 76 cols[0] = (i-1+n)%n; 77 cols[1] = i%n; 78 cols[2] = (i+1)%n; 79 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 80 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 81 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 82 ierr = MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES);CHKERRQ(ierr); 83 } 84 if (symmetric) { 85 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 86 } 87 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 88 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 89 90 ierr = MatISGetLocalMat(A,&B);CHKERRQ(ierr); 91 ierr = MatGetType(B,&lmtype);CHKERRQ(ierr); 92 93 /* test MatGetInfo */ 94 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n");CHKERRQ(ierr); 95 if (!PetscGlobalRank) { 96 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 97 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 98 } 99 ierr = MatGetInfo(A,MAT_LOCAL,&info);CHKERRQ(ierr); 100 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 101 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process %2d: %D %D %D %D %D\n",PetscGlobalRank,(PetscInt)info.nz_used, 102 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);CHKERRQ(ierr); 103 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 104 ierr = MatGetInfo(A,MAT_GLOBAL_MAX,&info);CHKERRQ(ierr); 105 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax : %D %D %D %D %D\n",(PetscInt)info.nz_used, 106 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);CHKERRQ(ierr); 107 ierr = MatGetInfo(A,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); 108 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum : %D %D %D %D %D\n",(PetscInt)info.nz_used, 109 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);CHKERRQ(ierr); 110 111 /* test MatView */ 112 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n");CHKERRQ(ierr); 113 ierr = MatView(A,NULL);CHKERRQ(ierr); 114 115 /* Create a MPIAIJ matrix, same as A */ 116 ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr); 117 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 118 ierr = MatSetType(B,MATAIJ);CHKERRQ(ierr); 119 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 120 ierr = MatSetUp(B);CHKERRQ(ierr); 121 ierr = MatSetLocalToGlobalMapping(B,rmap,cmap);CHKERRQ(ierr); 122 ierr = MatMPIAIJSetPreallocation(B,3,NULL,3,NULL);CHKERRQ(ierr); 123 ierr = MatMPIBAIJSetPreallocation(B,1,3,NULL,3,NULL);CHKERRQ(ierr); 124 #if defined(PETSC_HAVE_HYPRE) 125 ierr = MatHYPRESetPreallocation(B,3,NULL,3,NULL);CHKERRQ(ierr); 126 #endif 127 ierr = MatISSetPreallocation(B,3,NULL,3,NULL);CHKERRQ(ierr); 128 for (i=0; i<m; i++) { 129 PetscScalar v[3]; 130 PetscInt cols[3]; 131 132 cols[0] = (i-1+n)%n; 133 cols[1] = i%n; 134 cols[2] = (i+1)%n; 135 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 136 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 137 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 138 ierr = MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES);CHKERRQ(ierr); 139 } 140 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 143 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 144 145 /* test CheckMat */ 146 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test CheckMat\n");CHKERRQ(ierr); 147 ierr = CheckMat(A,B,PETSC_FALSE,"CheckMat");CHKERRQ(ierr); 148 149 /* test MatDuplicate and MatAXPY */ 150 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatDuplicate and MatAXPY\n");CHKERRQ(ierr); 151 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 152 ierr = CheckMat(A,A2,PETSC_FALSE,"MatDuplicate and MatAXPY");CHKERRQ(ierr); 153 154 /* test MatConvert */ 155 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_IS_XAIJ\n");CHKERRQ(ierr); 156 ierr = MatConvert(A2,MATAIJ,MAT_INITIAL_MATRIX,&B2);CHKERRQ(ierr); 157 ierr = CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INITIAL_MATRIX");CHKERRQ(ierr); 158 ierr = MatConvert(A2,MATAIJ,MAT_REUSE_MATRIX,&B2);CHKERRQ(ierr); 159 ierr = CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_REUSE_MATRIX");CHKERRQ(ierr); 160 ierr = MatConvert(A2,MATAIJ,MAT_INPLACE_MATRIX,&A2);CHKERRQ(ierr); 161 ierr = CheckMat(B,A2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INPLACE_MATRIX");CHKERRQ(ierr); 162 ierr = MatDestroy(&A2);CHKERRQ(ierr); 163 ierr = MatDestroy(&B2);CHKERRQ(ierr); 164 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_XAIJ_IS\n");CHKERRQ(ierr); 165 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 166 ierr = MatConvert(B2,MATIS,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr); 167 ierr = CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INITIAL_MATRIX");CHKERRQ(ierr); 168 ierr = MatConvert(B2,MATIS,MAT_REUSE_MATRIX,&A2);CHKERRQ(ierr); 169 ierr = CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_REUSE_MATRIX");CHKERRQ(ierr); 170 ierr = MatConvert(B2,MATIS,MAT_INPLACE_MATRIX,&B2);CHKERRQ(ierr); 171 ierr = CheckMat(A,B2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INPLACE_MATRIX");CHKERRQ(ierr); 172 ierr = MatDestroy(&A2);CHKERRQ(ierr); 173 ierr = MatDestroy(&B2);CHKERRQ(ierr); 174 ierr = PetscStrcmp(lmtype,MATSEQAIJ,&isaij);CHKERRQ(ierr); 175 if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */ 176 PetscInt ri, ci, rr[3] = {0,1,0}, cr[4] = {1,2,0,1}, rk[3] = {0,2,1}, ck[4] = {1,0,3,2}; 177 178 for (ri = 0; ri < 2; ri++) { 179 PetscInt *r; 180 181 r = (PetscInt*)(ri == 0 ? rr : rk); 182 for (ci = 0; ci < 2; ci++) { 183 PetscInt *c,rb,cb; 184 185 c = (PetscInt*)(ci == 0 ? cr : ck); 186 for (rb = 1; rb < 4; rb++) { 187 ierr = ISCreateBlock(PETSC_COMM_SELF,rb,3,r,PETSC_COPY_VALUES,&is);CHKERRQ(ierr); 188 ierr = ISLocalToGlobalMappingCreateIS(is,&rmap);CHKERRQ(ierr); 189 ierr = ISDestroy(&is);CHKERRQ(ierr); 190 for (cb = 1; cb < 4; cb++) { 191 Mat T,lT,T2; 192 char testname[256]; 193 194 ierr = PetscSNPrintf(testname,sizeof(testname),"MatConvert_IS_XAIJ special case (%D %D, bs %D %D)",ri,ci,rb,cb);CHKERRQ(ierr); 195 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test %s\n",testname);CHKERRQ(ierr); 196 197 ierr = ISCreateBlock(PETSC_COMM_SELF,cb,4,c,PETSC_COPY_VALUES,&is);CHKERRQ(ierr); 198 ierr = ISLocalToGlobalMappingCreateIS(is,&cmap);CHKERRQ(ierr); 199 ierr = ISDestroy(&is);CHKERRQ(ierr); 200 201 ierr = MatCreate(PETSC_COMM_SELF,&T);CHKERRQ(ierr); 202 ierr = MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,rb*3,cb*4);CHKERRQ(ierr); 203 ierr = MatSetType(T,MATIS);CHKERRQ(ierr); 204 ierr = MatSetLocalToGlobalMapping(T,rmap,cmap);CHKERRQ(ierr); 205 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 206 ierr = MatISGetLocalMat(T,&lT);CHKERRQ(ierr); 207 ierr = MatSetType(lT,MATSEQAIJ);CHKERRQ(ierr); 208 ierr = MatSeqAIJSetPreallocation(lT,cb*4,NULL);CHKERRQ(ierr); 209 ierr = MatSetRandom(lT,NULL);CHKERRQ(ierr); 210 ierr = MatConvert(lT,lmtype,MAT_INPLACE_MATRIX,&lT);CHKERRQ(ierr); 211 ierr = MatISRestoreLocalMat(T,&lT);CHKERRQ(ierr); 212 ierr = MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 ierr = MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 214 215 ierr = MatConvert(T,MATAIJ,MAT_INITIAL_MATRIX,&T2);CHKERRQ(ierr); 216 ierr = CheckMat(T,T2,PETSC_TRUE,"MAT_INITIAL_MATRIX");CHKERRQ(ierr); 217 ierr = MatConvert(T,MATAIJ,MAT_REUSE_MATRIX,&T2);CHKERRQ(ierr); 218 ierr = CheckMat(T,T2,PETSC_TRUE,"MAT_REUSE_MATRIX");CHKERRQ(ierr); 219 ierr = MatDestroy(&T2);CHKERRQ(ierr); 220 ierr = MatDuplicate(T,MAT_COPY_VALUES,&T2);CHKERRQ(ierr); 221 ierr = MatConvert(T2,MATAIJ,MAT_INPLACE_MATRIX,&T2);CHKERRQ(ierr); 222 ierr = CheckMat(T,T2,PETSC_TRUE,"MAT_INPLACE_MATRIX");CHKERRQ(ierr); 223 ierr = MatDestroy(&T);CHKERRQ(ierr); 224 ierr = MatDestroy(&T2);CHKERRQ(ierr); 225 } 226 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 227 } 228 } 229 } 230 } 231 232 /* test MatDiagonalScale */ 233 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n");CHKERRQ(ierr); 234 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 235 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 236 ierr = MatCreateVecs(A,&x,&y);CHKERRQ(ierr); 237 ierr = VecSetRandom(x,NULL);CHKERRQ(ierr); 238 if (symmetric) { 239 ierr = VecCopy(x,y);CHKERRQ(ierr); 240 } else { 241 ierr = VecSetRandom(y,NULL);CHKERRQ(ierr); 242 ierr = VecScale(y,8.);CHKERRQ(ierr); 243 } 244 ierr = MatDiagonalScale(A2,y,x);CHKERRQ(ierr); 245 ierr = MatDiagonalScale(B2,y,x);CHKERRQ(ierr); 246 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale");CHKERRQ(ierr); 247 ierr = MatDestroy(&A2);CHKERRQ(ierr); 248 ierr = MatDestroy(&B2);CHKERRQ(ierr); 249 ierr = VecDestroy(&x);CHKERRQ(ierr); 250 ierr = VecDestroy(&y);CHKERRQ(ierr); 251 252 /* test MatPtAP (A IS and B AIJ) */ 253 if (isaij && m == n) { 254 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatPtAP\n");CHKERRQ(ierr); 255 ierr = MatISStoreL2L(A,PETSC_TRUE);CHKERRQ(ierr); 256 ierr = MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&A2);CHKERRQ(ierr); 257 ierr = MatPtAP(B,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B2);CHKERRQ(ierr); 258 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_INITIAL_MATRIX");CHKERRQ(ierr); 259 ierr = MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&A2);CHKERRQ(ierr); 260 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_REUSE_MATRIX");CHKERRQ(ierr); 261 ierr = MatDestroy(&A2);CHKERRQ(ierr); 262 ierr = MatDestroy(&B2);CHKERRQ(ierr); 263 } 264 265 /* test MatGetLocalSubMatrix */ 266 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n");CHKERRQ(ierr); 267 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2);CHKERRQ(ierr); 268 ierr = ISCreateStride(PETSC_COMM_WORLD,m/2+m%2,0,2,&reven);CHKERRQ(ierr); 269 ierr = ISCreateStride(PETSC_COMM_WORLD,m/2,1,2,&rodd);CHKERRQ(ierr); 270 ierr = ISCreateStride(PETSC_COMM_WORLD,n/2+n%2,0,2,&ceven);CHKERRQ(ierr); 271 ierr = ISCreateStride(PETSC_COMM_WORLD,n/2,1,2,&codd);CHKERRQ(ierr); 272 ierr = MatGetLocalSubMatrix(A2,reven,ceven,&Aee);CHKERRQ(ierr); 273 ierr = MatGetLocalSubMatrix(A2,reven,codd,&Aeo);CHKERRQ(ierr); 274 ierr = MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe);CHKERRQ(ierr); 275 ierr = MatGetLocalSubMatrix(A2,rodd,codd,&Aoo);CHKERRQ(ierr); 276 for (i=0; i<m; i++) { 277 PetscInt j,je,jo,colse[3], colso[3]; 278 PetscScalar ve[3], vo[3]; 279 PetscScalar v[3]; 280 PetscInt cols[3]; 281 282 cols[0] = (i-1+n)%n; 283 cols[1] = i%n; 284 cols[2] = (i+1)%n; 285 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 286 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 287 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 288 for (j=0,je=0,jo=0;j<3;j++) { 289 if (cols[j]%2) { 290 vo[jo] = v[j]; 291 colso[jo++] = cols[j]/2; 292 } else { 293 ve[je] = v[j]; 294 colse[je++] = cols[j]/2; 295 } 296 } 297 if (i%2) { 298 PetscInt row = i/2; 299 ierr = MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES);CHKERRQ(ierr); 300 ierr = MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES);CHKERRQ(ierr); 301 } else { 302 PetscInt row = i/2; 303 ierr = MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES);CHKERRQ(ierr); 304 ierr = MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES);CHKERRQ(ierr); 305 } 306 } 307 ierr = MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee);CHKERRQ(ierr); 308 ierr = MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo);CHKERRQ(ierr); 309 ierr = MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe);CHKERRQ(ierr); 310 ierr = MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo);CHKERRQ(ierr); 311 ierr = ISDestroy(&reven);CHKERRQ(ierr); 312 ierr = ISDestroy(&ceven);CHKERRQ(ierr); 313 ierr = ISDestroy(&rodd);CHKERRQ(ierr); 314 ierr = ISDestroy(&codd);CHKERRQ(ierr); 315 ierr = MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 316 ierr = MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 317 ierr = MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 318 ierr = CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix");CHKERRQ(ierr); 319 ierr = MatDestroy(&A2);CHKERRQ(ierr); 320 321 /* test MatConvert_Nest_IS */ 322 testT = PETSC_FALSE; 323 ierr = PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL);CHKERRQ(ierr); 324 325 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n");CHKERRQ(ierr); 326 nr = 2; 327 nc = 2; 328 ierr = PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL);CHKERRQ(ierr); 329 ierr = PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);CHKERRQ(ierr); 330 if (testT) { 331 ierr = MatGetOwnershipRange(A,&cst,&cen);CHKERRQ(ierr); 332 ierr = MatGetOwnershipRangeColumn(A,&rst,&ren);CHKERRQ(ierr); 333 } else { 334 ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr); 335 ierr = MatGetOwnershipRangeColumn(A,&cst,&cen);CHKERRQ(ierr); 336 } 337 ierr = PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats);CHKERRQ(ierr); 338 for (i=0;i<nr*nc;i++) { 339 if (testT) { 340 ierr = MatCreateTranspose(A,&mats[i]);CHKERRQ(ierr); 341 ierr = MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc]);CHKERRQ(ierr); 342 } else { 343 ierr = MatDuplicate(A,MAT_COPY_VALUES,&mats[i]);CHKERRQ(ierr); 344 ierr = MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc]);CHKERRQ(ierr); 345 } 346 } 347 for (i=0;i<nr;i++) { 348 ierr = ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i]);CHKERRQ(ierr); 349 } 350 for (i=0;i<nc;i++) { 351 ierr = ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i]);CHKERRQ(ierr); 352 } 353 ierr = MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2);CHKERRQ(ierr); 354 ierr = MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2);CHKERRQ(ierr); 355 for (i=0;i<nr;i++) { 356 ierr = ISDestroy(&rows[i]);CHKERRQ(ierr); 357 } 358 for (i=0;i<nc;i++) { 359 ierr = ISDestroy(&cols[i]);CHKERRQ(ierr); 360 } 361 for (i=0;i<2*nr*nc;i++) { 362 ierr = MatDestroy(&mats[i]);CHKERRQ(ierr); 363 } 364 ierr = PetscFree3(rows,cols,mats);CHKERRQ(ierr); 365 ierr = MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T);CHKERRQ(ierr); 366 ierr = MatDestroy(&B2);CHKERRQ(ierr); 367 ierr = MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2);CHKERRQ(ierr); 368 ierr = CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INITIAL_MATRIX");CHKERRQ(ierr); 369 ierr = MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2);CHKERRQ(ierr); 370 ierr = CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_REUSE_MATRIX");CHKERRQ(ierr); 371 ierr = MatDestroy(&B2);CHKERRQ(ierr); 372 ierr = MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2);CHKERRQ(ierr); 373 ierr = CheckMat(A2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INPLACE_MATRIX");CHKERRQ(ierr); 374 ierr = MatDestroy(&T);CHKERRQ(ierr); 375 ierr = MatDestroy(&A2);CHKERRQ(ierr); 376 377 /* test MatCreateSubMatrix */ 378 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n");CHKERRQ(ierr); 379 if (!rank) { 380 ierr = ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is);CHKERRQ(ierr); 381 ierr = ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2);CHKERRQ(ierr); 382 } else if (rank == 1) { 383 ierr = ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);CHKERRQ(ierr); 384 if (n > 3) { 385 ierr = ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2);CHKERRQ(ierr); 386 } else { 387 ierr = ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);CHKERRQ(ierr); 388 } 389 } else if (rank == 2 && n > 4) { 390 ierr = ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);CHKERRQ(ierr); 391 ierr = ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2);CHKERRQ(ierr); 392 } else { 393 ierr = ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);CHKERRQ(ierr); 394 ierr = ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);CHKERRQ(ierr); 395 } 396 ierr = MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr); 397 ierr = MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2);CHKERRQ(ierr); 398 ierr = CheckMat(A2,B2,PETSC_TRUE,"first MatCreateSubMatrix");CHKERRQ(ierr); 399 400 ierr = MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2);CHKERRQ(ierr); 401 ierr = MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2);CHKERRQ(ierr); 402 ierr = CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix");CHKERRQ(ierr); 403 ierr = MatDestroy(&A2);CHKERRQ(ierr); 404 ierr = MatDestroy(&B2);CHKERRQ(ierr); 405 406 if (!symmetric) { 407 ierr = MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr); 408 ierr = MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2);CHKERRQ(ierr); 409 ierr = MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2);CHKERRQ(ierr); 410 ierr = MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2);CHKERRQ(ierr); 411 ierr = CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix");CHKERRQ(ierr); 412 } 413 414 ierr = MatDestroy(&A2);CHKERRQ(ierr); 415 ierr = MatDestroy(&B2);CHKERRQ(ierr); 416 ierr = ISDestroy(&is);CHKERRQ(ierr); 417 ierr = ISDestroy(&is2);CHKERRQ(ierr); 418 419 /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */ 420 if (size > 1) { 421 if (!rank) { 422 PetscInt st,len; 423 424 st = (m+1)/2; 425 len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0)); 426 ierr = ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is);CHKERRQ(ierr); 427 } else { 428 ierr = ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);CHKERRQ(ierr); 429 } 430 } else { 431 ierr = ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);CHKERRQ(ierr); 432 } 433 434 if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */ 435 /* test MatDiagonalSet */ 436 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n");CHKERRQ(ierr); 437 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 438 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 439 ierr = MatCreateVecs(A,NULL,&x);CHKERRQ(ierr); 440 ierr = VecSetRandom(x,NULL);CHKERRQ(ierr); 441 ierr = MatDiagonalSet(A2,x,INSERT_VALUES);CHKERRQ(ierr); 442 ierr = MatDiagonalSet(B2,x,INSERT_VALUES);CHKERRQ(ierr); 443 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet");CHKERRQ(ierr); 444 ierr = VecDestroy(&x);CHKERRQ(ierr); 445 ierr = MatDestroy(&A2);CHKERRQ(ierr); 446 ierr = MatDestroy(&B2);CHKERRQ(ierr); 447 448 /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */ 449 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n");CHKERRQ(ierr); 450 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 451 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 452 ierr = MatShift(A2,2.0);CHKERRQ(ierr); 453 ierr = MatShift(B2,2.0);CHKERRQ(ierr); 454 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatShift");CHKERRQ(ierr); 455 ierr = MatDestroy(&A2);CHKERRQ(ierr); 456 ierr = MatDestroy(&B2);CHKERRQ(ierr); 457 458 /* nonzero diag value is supported for square matrices only */ 459 ierr = TestMatZeroRows(A,B,PETSC_TRUE,is,diag);CHKERRQ(ierr); 460 } 461 ierr = TestMatZeroRows(A,B,squaretest,is,0.0);CHKERRQ(ierr); 462 ierr = ISDestroy(&is);CHKERRQ(ierr); 463 464 /* test MatTranspose */ 465 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n");CHKERRQ(ierr); 466 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr); 467 ierr = MatTranspose(B,MAT_INITIAL_MATRIX,&B2);CHKERRQ(ierr); 468 ierr = CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose");CHKERRQ(ierr); 469 470 ierr = MatTranspose(A,MAT_REUSE_MATRIX,&A2);CHKERRQ(ierr); 471 ierr = CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose");CHKERRQ(ierr); 472 ierr = MatDestroy(&A2);CHKERRQ(ierr); 473 474 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 475 ierr = MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);CHKERRQ(ierr); 476 ierr = CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose");CHKERRQ(ierr); 477 ierr = MatDestroy(&A2);CHKERRQ(ierr); 478 479 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr); 480 ierr = CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose");CHKERRQ(ierr); 481 ierr = MatDestroy(&A2);CHKERRQ(ierr); 482 ierr = MatDestroy(&B2);CHKERRQ(ierr); 483 484 /* test MatISFixLocalEmpty */ 485 if (isaij) { 486 PetscInt r[2]; 487 488 r[0] = 0; 489 r[1] = PetscMin(m,n)-1; 490 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatISFixLocalEmpty\n");CHKERRQ(ierr); 491 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 492 ierr = MatISFixLocalEmpty(A2,PETSC_TRUE);CHKERRQ(ierr); 493 ierr = MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 494 ierr = MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 495 ierr = CheckMat(A2,B,PETSC_FALSE,"MatISFixLocalEmpty (null)");CHKERRQ(ierr); 496 497 ierr = MatZeroRows(A2,2,r,0.0,NULL,NULL);CHKERRQ(ierr); 498 ierr = MatViewFromOptions(A2,NULL,"-fixempty_view");CHKERRQ(ierr); 499 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 500 ierr = MatZeroRows(B2,2,r,0.0,NULL,NULL);CHKERRQ(ierr); 501 ierr = MatISFixLocalEmpty(A2,PETSC_TRUE);CHKERRQ(ierr); 502 ierr = MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 503 ierr = MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 504 ierr = MatViewFromOptions(A2,NULL,"-fixempty_view");CHKERRQ(ierr); 505 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows)");CHKERRQ(ierr); 506 ierr = MatDestroy(&A2);CHKERRQ(ierr); 507 508 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 509 ierr = MatZeroRows(A2,2,r,0.0,NULL,NULL);CHKERRQ(ierr); 510 ierr = MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);CHKERRQ(ierr); 511 ierr = MatTranspose(B2,MAT_INPLACE_MATRIX,&B2);CHKERRQ(ierr); 512 ierr = MatViewFromOptions(A2,NULL,"-fixempty_view");CHKERRQ(ierr); 513 ierr = MatISFixLocalEmpty(A2,PETSC_TRUE);CHKERRQ(ierr); 514 ierr = MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 515 ierr = MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 516 ierr = MatViewFromOptions(A2,NULL,"-fixempty_view");CHKERRQ(ierr); 517 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (cols)");CHKERRQ(ierr); 518 519 ierr = MatDestroy(&A2);CHKERRQ(ierr); 520 ierr = MatDestroy(&B2);CHKERRQ(ierr); 521 522 if (squaretest) { 523 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 524 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 525 ierr = MatZeroRowsColumns(A2,2,r,0.0,NULL,NULL);CHKERRQ(ierr); 526 ierr = MatZeroRowsColumns(B2,2,r,0.0,NULL,NULL);CHKERRQ(ierr); 527 ierr = MatViewFromOptions(A2,NULL,"-fixempty_view");CHKERRQ(ierr); 528 ierr = MatISFixLocalEmpty(A2,PETSC_TRUE);CHKERRQ(ierr); 529 ierr = MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 530 ierr = MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 531 ierr = MatViewFromOptions(A2,NULL,"-fixempty_view");CHKERRQ(ierr); 532 ierr = CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows+cols)");CHKERRQ(ierr); 533 ierr = MatDestroy(&A2);CHKERRQ(ierr); 534 ierr = MatDestroy(&B2);CHKERRQ(ierr); 535 } 536 537 } 538 539 /* test MatInvertBlockDiagonal 540 special cases for block-diagonal matrices */ 541 if (m == n) { 542 ISLocalToGlobalMapping map; 543 Mat Abd,Bbd; 544 IS is,bis; 545 const PetscScalar *isbd,*aijbd; 546 PetscScalar *vals; 547 const PetscInt *sts,*idxs; 548 PetscInt *idxs2,diff,perm,nl,bs,st,en,in; 549 PetscBool ok; 550 551 for (diff = 0; diff < 3; diff++) { 552 for (perm = 0; perm < 3; perm++) { 553 for (bs = 1; bs < 4; bs++) { 554 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatInvertBlockDiagonal blockdiag %D %D %D %D\n",n,diff,perm,bs);CHKERRQ(ierr); 555 ierr = PetscMalloc1(bs*bs,&vals);CHKERRQ(ierr); 556 ierr = MatGetOwnershipRanges(A,&sts);CHKERRQ(ierr); 557 switch (diff) { 558 case 1: /* inverted layout by processes */ 559 in = 1; 560 st = sts[size - rank - 1]; 561 en = sts[size - rank]; 562 nl = en - st; 563 break; 564 case 2: /* round-robin layout */ 565 in = size; 566 st = rank; 567 nl = n/size; 568 if (rank < n%size) nl++; 569 break; 570 default: /* same layout */ 571 in = 1; 572 st = sts[rank]; 573 en = sts[rank + 1]; 574 nl = en - st; 575 break; 576 } 577 ierr = ISCreateStride(PETSC_COMM_WORLD,nl,st,in,&is);CHKERRQ(ierr); 578 ierr = ISGetLocalSize(is,&nl);CHKERRQ(ierr); 579 ierr = ISGetIndices(is,&idxs);CHKERRQ(ierr); 580 ierr = PetscMalloc1(nl,&idxs2);CHKERRQ(ierr); 581 for (i=0;i<nl;i++) { 582 switch (perm) { /* invert some of the indices */ 583 case 2: 584 idxs2[i] = rank%2 ? idxs[i] : idxs[nl-i-1]; 585 break; 586 case 1: 587 idxs2[i] = rank%2 ? idxs[nl-i-1] : idxs[i]; 588 break; 589 default: 590 idxs2[i] = idxs[i]; 591 break; 592 } 593 } 594 ierr = ISRestoreIndices(is,&idxs);CHKERRQ(ierr); 595 ierr = ISCreateBlock(PETSC_COMM_WORLD,bs,nl,idxs2,PETSC_OWN_POINTER,&bis);CHKERRQ(ierr); 596 ierr = ISLocalToGlobalMappingCreateIS(bis,&map);CHKERRQ(ierr); 597 ierr = MatCreateIS(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,bs*n,bs*n,map,NULL,&Abd);CHKERRQ(ierr); 598 ierr = ISLocalToGlobalMappingDestroy(&map);CHKERRQ(ierr); 599 ierr = MatISSetPreallocation(Abd,bs,NULL,0,NULL);CHKERRQ(ierr); 600 for (i=0;i<nl;i++) { 601 PetscInt b1,b2; 602 603 for (b1=0;b1<bs;b1++) { 604 for (b2=0;b2<bs;b2++) { 605 vals[b1*bs + b2] = i*bs*bs + b1*bs + b2 + 1 + (b1 == b2 ? 1.0 : 0); 606 } 607 } 608 ierr = MatSetValuesBlockedLocal(Abd,1,&i,1,&i,vals,INSERT_VALUES);CHKERRQ(ierr); 609 } 610 ierr = MatAssemblyBegin(Abd,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 611 ierr = MatAssemblyEnd(Abd,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 612 ierr = MatConvert(Abd,MATAIJ,MAT_INITIAL_MATRIX,&Bbd);CHKERRQ(ierr); 613 ierr = MatInvertBlockDiagonal(Abd,&isbd);CHKERRQ(ierr); 614 ierr = MatInvertBlockDiagonal(Bbd,&aijbd);CHKERRQ(ierr); 615 ierr = MatGetLocalSize(Bbd,&nl,NULL);CHKERRQ(ierr); 616 ok = PETSC_TRUE; 617 for (i=0;i<nl/bs;i++) { 618 PetscInt b1,b2; 619 620 for (b1=0;b1<bs;b1++) { 621 for (b2=0;b2<bs;b2++) { 622 if (PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]-aijbd[i*bs*bs+b1*bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE; 623 if (!ok) { 624 ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] ERROR block %d, entry %d %d: %g %g\n",rank,i,b1,b2,isbd[i*bs*bs+b1*bs + b2],aijbd[i*bs*bs+b1*bs + b2]);CHKERRQ(ierr); 625 break; 626 } 627 } 628 if (!ok) break; 629 } 630 if (!ok) break; 631 } 632 ierr = MatDestroy(&Abd);CHKERRQ(ierr); 633 ierr = MatDestroy(&Bbd);CHKERRQ(ierr); 634 ierr = PetscFree(vals);CHKERRQ(ierr); 635 ierr = ISDestroy(&is);CHKERRQ(ierr); 636 ierr = ISDestroy(&bis);CHKERRQ(ierr); 637 } 638 } 639 } 640 } 641 /* free testing matrices */ 642 ierr = MatDestroy(&A);CHKERRQ(ierr); 643 ierr = MatDestroy(&B);CHKERRQ(ierr); 644 ierr = PetscFinalize(); 645 return ierr; 646 } 647 648 PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func) 649 { 650 Mat Bcheck; 651 PetscReal error; 652 PetscErrorCode ierr; 653 654 PetscFunctionBeginUser; 655 if (!usemult && B) { 656 PetscBool hasnorm; 657 658 ierr = MatHasOperation(B,MATOP_NORM,&hasnorm);CHKERRQ(ierr); 659 if (!hasnorm) usemult = PETSC_TRUE; 660 } 661 if (!usemult) { 662 if (B) { 663 MatType Btype; 664 665 ierr = MatGetType(B,&Btype);CHKERRQ(ierr); 666 ierr = MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck);CHKERRQ(ierr); 667 } else { 668 ierr = MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);CHKERRQ(ierr); 669 } 670 if (B) { /* if B is present, subtract it */ 671 ierr = MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 672 } 673 ierr = MatNorm(Bcheck,NORM_INFINITY,&error);CHKERRQ(ierr); 674 if (error > PETSC_SQRT_MACHINE_EPSILON) { 675 ISLocalToGlobalMapping rl2g,cl2g; 676 677 ierr = PetscObjectSetName((PetscObject)Bcheck,"Assembled Bcheck");CHKERRQ(ierr); 678 ierr = MatView(Bcheck,NULL);CHKERRQ(ierr); 679 if (B) { 680 ierr = PetscObjectSetName((PetscObject)B,"Assembled AIJ");CHKERRQ(ierr); 681 ierr = MatView(B,NULL);CHKERRQ(ierr); 682 ierr = MatDestroy(&Bcheck);CHKERRQ(ierr); 683 ierr = MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);CHKERRQ(ierr); 684 ierr = PetscObjectSetName((PetscObject)Bcheck,"Assembled IS");CHKERRQ(ierr); 685 ierr = MatView(Bcheck,NULL);CHKERRQ(ierr); 686 } 687 ierr = MatDestroy(&Bcheck);CHKERRQ(ierr); 688 ierr = PetscObjectSetName((PetscObject)A,"MatIS");CHKERRQ(ierr); 689 ierr = MatView(A,NULL);CHKERRQ(ierr); 690 ierr = MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);CHKERRQ(ierr); 691 ierr = ISLocalToGlobalMappingView(rl2g,NULL);CHKERRQ(ierr); 692 ierr = ISLocalToGlobalMappingView(cl2g,NULL);CHKERRQ(ierr); 693 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,error); 694 } 695 ierr = MatDestroy(&Bcheck);CHKERRQ(ierr); 696 } else { 697 PetscBool ok,okt; 698 699 ierr = MatMultEqual(A,B,3,&ok);CHKERRQ(ierr); 700 ierr = MatMultTransposeEqual(A,B,3,&okt);CHKERRQ(ierr); 701 if (!ok || !okt) SETERRQ3(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: mult ok ? %d, multtranspose ok ? %d",func,ok,okt); 702 } 703 PetscFunctionReturn(0); 704 } 705 706 PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag) 707 { 708 Mat B,Bcheck,B2 = NULL,lB; 709 Vec x = NULL, b = NULL, b2 = NULL; 710 ISLocalToGlobalMapping l2gr,l2gc; 711 PetscReal error; 712 char diagstr[16]; 713 const PetscInt *idxs; 714 PetscInt rst,ren,i,n,N,d; 715 PetscMPIInt rank; 716 PetscBool miss,haszerorows; 717 PetscErrorCode ierr; 718 719 PetscFunctionBeginUser; 720 if (diag == 0.) { 721 ierr = PetscStrcpy(diagstr,"zero");CHKERRQ(ierr); 722 } else { 723 ierr = PetscStrcpy(diagstr,"nonzero");CHKERRQ(ierr); 724 } 725 ierr = ISView(is,NULL);CHKERRQ(ierr); 726 ierr = MatGetLocalToGlobalMapping(A,&l2gr,&l2gc);CHKERRQ(ierr); 727 /* tests MatDuplicate and MatCopy */ 728 if (diag == 0.) { 729 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 730 } else { 731 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);CHKERRQ(ierr); 732 ierr = MatCopy(A,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 733 } 734 ierr = MatISGetLocalMat(B,&lB);CHKERRQ(ierr); 735 ierr = MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows);CHKERRQ(ierr); 736 if (squaretest && haszerorows) { 737 738 ierr = MatCreateVecs(B,&x,&b);CHKERRQ(ierr); 739 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 740 ierr = VecSetLocalToGlobalMapping(b,l2gr);CHKERRQ(ierr); 741 ierr = VecSetLocalToGlobalMapping(x,l2gc);CHKERRQ(ierr); 742 ierr = VecSetRandom(x,NULL);CHKERRQ(ierr); 743 ierr = VecSetRandom(b,NULL);CHKERRQ(ierr); 744 /* mimic b[is] = x[is] */ 745 ierr = VecDuplicate(b,&b2);CHKERRQ(ierr); 746 ierr = VecSetLocalToGlobalMapping(b2,l2gr);CHKERRQ(ierr); 747 ierr = VecCopy(b,b2);CHKERRQ(ierr); 748 ierr = ISGetLocalSize(is,&n);CHKERRQ(ierr); 749 ierr = ISGetIndices(is,&idxs);CHKERRQ(ierr); 750 ierr = VecGetSize(x,&N);CHKERRQ(ierr); 751 for (i=0;i<n;i++) { 752 if (0 <= idxs[i] && idxs[i] < N) { 753 ierr = VecSetValue(b2,idxs[i],diag,INSERT_VALUES);CHKERRQ(ierr); 754 ierr = VecSetValue(x,idxs[i],1.,INSERT_VALUES);CHKERRQ(ierr); 755 } 756 } 757 ierr = VecAssemblyBegin(b2);CHKERRQ(ierr); 758 ierr = VecAssemblyEnd(b2);CHKERRQ(ierr); 759 ierr = VecAssemblyBegin(x);CHKERRQ(ierr); 760 ierr = VecAssemblyEnd(x);CHKERRQ(ierr); 761 ierr = ISRestoreIndices(is,&idxs);CHKERRQ(ierr); 762 /* test ZeroRows on MATIS */ 763 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);CHKERRQ(ierr); 764 ierr = MatZeroRowsIS(B,is,diag,x,b);CHKERRQ(ierr); 765 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr);CHKERRQ(ierr); 766 ierr = MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL);CHKERRQ(ierr); 767 } else if (haszerorows) { 768 /* test ZeroRows on MATIS */ 769 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);CHKERRQ(ierr); 770 ierr = MatZeroRowsIS(B,is,diag,NULL,NULL);CHKERRQ(ierr); 771 b = b2 = x = NULL; 772 } else { 773 ierr = PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr);CHKERRQ(ierr); 774 b = b2 = x = NULL; 775 } 776 777 if (squaretest && haszerorows) { 778 ierr = VecAXPY(b2,-1.,b);CHKERRQ(ierr); 779 ierr = VecNorm(b2,NORM_INFINITY,&error);CHKERRQ(ierr); 780 if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",error,diagstr); 781 } 782 783 /* test MatMissingDiagonal */ 784 ierr = PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n");CHKERRQ(ierr); 785 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 786 ierr = MatMissingDiagonal(B,&miss,&d);CHKERRQ(ierr); 787 ierr = MatGetOwnershipRange(B,&rst,&ren);CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 789 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%D,%D) Missing %d, row %D (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr);CHKERRQ(ierr); 790 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 791 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 792 793 ierr = VecDestroy(&x);CHKERRQ(ierr); 794 ierr = VecDestroy(&b);CHKERRQ(ierr); 795 ierr = VecDestroy(&b2);CHKERRQ(ierr); 796 797 /* check the result of ZeroRows with that from MPIAIJ routines 798 assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */ 799 if (haszerorows) { 800 ierr = MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck);CHKERRQ(ierr); 801 ierr = MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 802 ierr = MatZeroRowsIS(Bcheck,is,diag,NULL,NULL);CHKERRQ(ierr); 803 ierr = CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows");CHKERRQ(ierr); 804 ierr = MatDestroy(&Bcheck);CHKERRQ(ierr); 805 } 806 ierr = MatDestroy(&B);CHKERRQ(ierr); 807 808 if (B2) { /* test MatZeroRowsColumns */ 809 ierr = MatDuplicate(Afull,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 810 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 811 ierr = MatZeroRowsColumnsIS(B,is,diag,NULL,NULL);CHKERRQ(ierr); 812 ierr = CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns");CHKERRQ(ierr); 813 ierr = MatDestroy(&B);CHKERRQ(ierr); 814 ierr = MatDestroy(&B2);CHKERRQ(ierr); 815 } 816 PetscFunctionReturn(0); 817 } 818 819 820 /*TEST 821 822 test: 823 args: -test_trans 824 825 test: 826 suffix: 2 827 nsize: 4 828 args: -matis_convert_local_nest -nr 3 -nc 4 829 830 test: 831 suffix: 3 832 nsize: 5 833 args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1 834 835 test: 836 suffix: 4 837 nsize: 6 838 args: -m 9 -n 12 -test_trans -nr 2 -nc 7 839 840 test: 841 suffix: 5 842 nsize: 6 843 args: -m 12 -n 12 -test_trans -nr 3 -nc 1 844 845 test: 846 suffix: 6 847 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 848 849 test: 850 suffix: 7 851 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 852 853 test: 854 suffix: 8 855 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 856 857 test: 858 suffix: 9 859 nsize: 5 860 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 861 862 test: 863 suffix: 10 864 nsize: 5 865 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 866 867 test: 868 suffix: vscat_default 869 nsize: 5 870 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 871 output_file: output/ex23_11.out 872 873 test: 874 suffix: 12 875 nsize: 3 876 args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3 877 878 testset: 879 output_file: output/ex23_13.out 880 nsize: 3 881 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap 882 filter: grep -v "type:" 883 test: 884 suffix: baij 885 args: -matis_localmat_type baij 886 test: 887 requires: viennacl 888 suffix: viennacl 889 args: -matis_localmat_type aijviennacl 890 test: 891 requires: cuda 892 suffix: cusparse 893 args: -matis_localmat_type aijcusparse 894 895 TEST*/ 896