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