1 2 static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Vec x,y,b,s1,s2; 9 Mat A; /* linear system matrix */ 10 Mat sA; /* symmetric part of the matrices */ 11 PetscInt n,mbs=16,bs=1,nz=3,prob=2,i,j,col[3],row,Ii,J,n1; 12 const PetscInt *ip_ptr; 13 PetscScalar neg_one = -1.0,value[3],alpha=0.1; 14 PetscMPIInt size; 15 IS ip, isrow, iscol; 16 PetscRandom rdm; 17 PetscBool reorder=PETSC_FALSE; 18 MatInfo minfo1,minfo2; 19 PetscReal norm1,norm2,tol=10*PETSC_SMALL; 20 21 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 22 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 23 PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 24 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL)); 25 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL)); 26 27 n = mbs*bs; 28 CHKERRQ(MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &A)); 29 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 30 CHKERRQ(MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &sA)); 31 CHKERRQ(MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 32 33 /* Test MatGetOwnershipRange() */ 34 CHKERRQ(MatGetOwnershipRange(A,&Ii,&J)); 35 CHKERRQ(MatGetOwnershipRange(sA,&i,&j)); 36 if (i-Ii || j-J) { 37 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n")); 38 } 39 40 /* Assemble matrix */ 41 if (bs == 1) { 42 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL)); 43 if (prob == 1) { /* tridiagonal matrix */ 44 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 45 for (i=1; i<n-1; i++) { 46 col[0] = i-1; col[1] = i; col[2] = i+1; 47 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 48 CHKERRQ(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 49 } 50 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1; 51 52 value[0]= 0.1; value[1]=-1; value[2]=2; 53 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 54 CHKERRQ(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 55 56 i = 0; col[0] = 0; col[1] = 1; col[2]=n-1; 57 58 value[0] = 2.0; value[1] = -1.0; value[2]=0.1; 59 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 60 CHKERRQ(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 61 } else if (prob ==2) { /* matrix for the five point stencil */ 62 n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001); 63 PetscCheckFalse(n1*n1 - n,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive integer!"); 64 for (i=0; i<n1; i++) { 65 for (j=0; j<n1; j++) { 66 Ii = j + n1*i; 67 if (i>0) { 68 J = Ii - n1; 69 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 70 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 71 } 72 if (i<n1-1) { 73 J = Ii + n1; 74 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 75 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 76 } 77 if (j>0) { 78 J = Ii - 1; 79 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 80 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 81 } 82 if (j<n1-1) { 83 J = Ii + 1; 84 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 85 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 86 } 87 } 88 } 89 } 90 } else { /* bs > 1 */ 91 #if defined(DIAGB) 92 for (block=0; block<n/bs; block++) { 93 /* diagonal blocks */ 94 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0; 95 for (i=1+block*bs; i<bs-1+block*bs; i++) { 96 col[0] = i-1; col[1] = i; col[2] = i+1; 97 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 98 CHKERRQ(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 99 } 100 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs; 101 102 value[0]=-1.0; value[1]=4.0; 103 CHKERRQ(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 104 CHKERRQ(MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES)); 105 106 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs; 107 108 value[0]=4.0; value[1] = -1.0; 109 CHKERRQ(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 110 CHKERRQ(MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES)); 111 } 112 #endif 113 /* off-diagonal blocks */ 114 value[0]=-1.0; 115 for (i=0; i<(n/bs-1)*bs; i++) { 116 col[0]=i+bs; 117 CHKERRQ(MatSetValues(A,1,&i,1,col,value,INSERT_VALUES)); 118 CHKERRQ(MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES)); 119 col[0]=i; row=i+bs; 120 CHKERRQ(MatSetValues(A,1,&row,1,col,value,INSERT_VALUES)); 121 CHKERRQ(MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES)); 122 } 123 } 124 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 125 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 126 127 CHKERRQ(MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY)); 128 CHKERRQ(MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY)); 129 130 /* Test MatNorm() */ 131 CHKERRQ(MatNorm(A,NORM_FROBENIUS,&norm1)); 132 CHKERRQ(MatNorm(sA,NORM_FROBENIUS,&norm2)); 133 norm1 -= norm2; 134 if (norm1<-tol || norm1>tol) { 135 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",(double)norm1)); 136 } 137 CHKERRQ(MatNorm(A,NORM_INFINITY,&norm1)); 138 CHKERRQ(MatNorm(sA,NORM_INFINITY,&norm2)); 139 norm1 -= norm2; 140 if (norm1<-tol || norm1>tol) { 141 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",(double)norm1)); 142 } 143 144 /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */ 145 CHKERRQ(MatGetInfo(A,MAT_LOCAL,&minfo1)); 146 CHKERRQ(MatGetInfo(sA,MAT_LOCAL,&minfo2)); 147 i = (int) (minfo1.nz_used - minfo2.nz_used); 148 j = (int) (minfo1.nz_allocated - minfo2.nz_allocated); 149 if (i<0 || j<0) { 150 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n")); 151 } 152 153 CHKERRQ(MatGetSize(A,&Ii,&J)); 154 CHKERRQ(MatGetSize(sA,&i,&j)); 155 if (i-Ii || j-J) { 156 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n")); 157 } 158 159 CHKERRQ(MatGetBlockSize(A, &Ii)); 160 CHKERRQ(MatGetBlockSize(sA, &i)); 161 if (i-Ii) { 162 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n")); 163 } 164 165 /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */ 166 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 167 CHKERRQ(PetscRandomSetFromOptions(rdm)); 168 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 169 CHKERRQ(VecDuplicate(x,&s1)); 170 CHKERRQ(VecDuplicate(x,&s2)); 171 CHKERRQ(VecDuplicate(x,&y)); 172 CHKERRQ(VecDuplicate(x,&b)); 173 174 CHKERRQ(VecSetRandom(x,rdm)); 175 176 CHKERRQ(MatDiagonalScale(A,x,x)); 177 CHKERRQ(MatDiagonalScale(sA,x,x)); 178 179 CHKERRQ(MatGetDiagonal(A,s1)); 180 CHKERRQ(MatGetDiagonal(sA,s2)); 181 CHKERRQ(VecNorm(s1,NORM_1,&norm1)); 182 CHKERRQ(VecNorm(s2,NORM_1,&norm2)); 183 norm1 -= norm2; 184 if (norm1<-tol || norm1>tol) { 185 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n")); 186 } 187 188 CHKERRQ(MatScale(A,alpha)); 189 CHKERRQ(MatScale(sA,alpha)); 190 191 /* Test MatMult(), MatMultAdd() */ 192 for (i=0; i<40; i++) { 193 CHKERRQ(VecSetRandom(x,rdm)); 194 CHKERRQ(MatMult(A,x,s1)); 195 CHKERRQ(MatMult(sA,x,s2)); 196 CHKERRQ(VecNorm(s1,NORM_1,&norm1)); 197 CHKERRQ(VecNorm(s2,NORM_1,&norm2)); 198 norm1 -= norm2; 199 if (norm1<-tol || norm1>tol) { 200 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n")); 201 } 202 } 203 204 for (i=0; i<40; i++) { 205 CHKERRQ(VecSetRandom(x,rdm)); 206 CHKERRQ(VecSetRandom(y,rdm)); 207 CHKERRQ(MatMultAdd(A,x,y,s1)); 208 CHKERRQ(MatMultAdd(sA,x,y,s2)); 209 CHKERRQ(VecNorm(s1,NORM_1,&norm1)); 210 CHKERRQ(VecNorm(s2,NORM_1,&norm2)); 211 norm1 -= norm2; 212 if (norm1<-tol || norm1>tol) { 213 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n")); 214 } 215 } 216 217 /* Test MatReordering() */ 218 CHKERRQ(MatGetOrdering(A,MATORDERINGNATURAL,&isrow,&iscol)); 219 ip = isrow; 220 221 if (reorder) { 222 IS nip; 223 PetscInt *nip_ptr; 224 CHKERRQ(PetscMalloc1(mbs,&nip_ptr)); 225 CHKERRQ(ISGetIndices(ip,&ip_ptr)); 226 CHKERRQ(PetscArraycpy(nip_ptr,ip_ptr,mbs)); 227 i = nip_ptr[1]; nip_ptr[1] = nip_ptr[mbs-2]; nip_ptr[mbs-2] = i; 228 i = nip_ptr[0]; nip_ptr[0] = nip_ptr[mbs-1]; nip_ptr[mbs-1] = i; 229 CHKERRQ(ISRestoreIndices(ip,&ip_ptr)); 230 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,mbs,nip_ptr,PETSC_COPY_VALUES,&nip)); 231 CHKERRQ(PetscFree(nip_ptr)); 232 233 CHKERRQ(MatReorderingSeqSBAIJ(sA, ip)); 234 CHKERRQ(ISDestroy(&nip)); 235 } 236 237 CHKERRQ(ISDestroy(&iscol)); 238 CHKERRQ(ISDestroy(&isrow)); 239 CHKERRQ(MatDestroy(&A)); 240 CHKERRQ(MatDestroy(&sA)); 241 CHKERRQ(VecDestroy(&x)); 242 CHKERRQ(VecDestroy(&y)); 243 CHKERRQ(VecDestroy(&s1)); 244 CHKERRQ(VecDestroy(&s2)); 245 CHKERRQ(VecDestroy(&b)); 246 CHKERRQ(PetscRandomDestroy(&rdm)); 247 248 CHKERRQ(PetscFinalize()); 249 return 0; 250 } 251 252 /*TEST 253 254 test: 255 args: -bs {{1 2 3 4 5 6 7 8}} 256 257 TEST*/ 258