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