1 2 static char help[] = "Tests MatInvertVariableBlockEnvelope()\n\n"; 3 4 #include <petscmat.h> 5 extern PetscErrorCode MatIsDiagonal(Mat); 6 extern PetscErrorCode BuildMatrix(const PetscInt*, PetscInt,const PetscInt *,Mat*); 7 8 int main(int argc,char **argv) 9 { 10 Mat A,C,D,F; 11 PetscInt i,j,rows[2],*parts,cnt, N = 21,nblocks,*blocksizes; 12 PetscScalar values[2][2]; 13 PetscReal rand; 14 PetscRandom rctx; 15 PetscMPIInt size; 16 17 PetscCall(PetscInitialize(&argc,&argv,(char*) 0,help)); 18 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE)); 19 20 PetscCall(MatCreate(PETSC_COMM_WORLD,&C)); 21 PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,6,18)); 22 PetscCall(MatSetFromOptions(C)); 23 PetscCall(MatSetUp(C)); 24 values[0][0] = 2; values[0][1] = 1; 25 values[1][0] = 1; values[1][1] = 2; 26 for (i=0;i<3;i++){ 27 rows[0] = 2*i; rows[1] = 2*i + 1; 28 PetscCall(MatSetValues(C,2,rows,2,rows,(PetscScalar*)values,INSERT_VALUES)); 29 } 30 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 31 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 32 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 33 34 PetscCall(MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DETERMINE,&A)); 35 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 36 37 PetscCall(MatInvertVariableBlockEnvelope(A,MAT_INITIAL_MATRIX,&D)); 38 PetscCall(MatView(D,PETSC_VIEWER_STDOUT_WORLD)); 39 40 PetscCall(MatMatMult(A,D,MAT_INITIAL_MATRIX,1.0,&F)); 41 PetscCall(MatView(F,PETSC_VIEWER_STDOUT_WORLD)); 42 PetscCall(MatIsDiagonal(F)); 43 44 PetscCall(MatDestroy(&A)); 45 PetscCall(MatDestroy(&D)); 46 PetscCall(MatDestroy(&C)); 47 PetscCall(MatDestroy(&F)); 48 49 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rctx)); 50 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 51 PetscCall(PetscMalloc1(size,&parts)); 52 53 for (j=0; j<128; j++) { 54 cnt = 0; 55 for (i=0; i<size-1; i++) { 56 PetscCall(PetscRandomGetValueReal(rctx,&rand)); 57 parts[i] = (PetscInt) N*rand; 58 parts[i] = PetscMin(parts[i],N-cnt); 59 cnt += parts[i]; 60 } 61 parts[size-1] = N - cnt; 62 63 PetscCall(PetscRandomGetValueReal(rctx,&rand)); 64 nblocks = rand*10; 65 nblocks = PetscMax(nblocks,2); 66 cnt = 0; 67 PetscCall(PetscMalloc1(nblocks,&blocksizes)); 68 for (i=0; i<nblocks-1; i++) { 69 PetscCall(PetscRandomGetValueReal(rctx,&rand)); 70 blocksizes[i] = PetscMax(1,(PetscInt) N*rand); 71 blocksizes[i] = PetscMin(blocksizes[i],N-cnt); 72 cnt += blocksizes[i]; 73 if (cnt == N) { 74 nblocks = i + 1; 75 break; 76 } 77 } 78 if (cnt < N) { 79 blocksizes[nblocks-1] = N - cnt; 80 } 81 82 PetscCall(BuildMatrix(parts,nblocks,blocksizes,&A)); 83 PetscCall(PetscFree(blocksizes)); 84 85 PetscCall(MatInvertVariableBlockEnvelope(A,MAT_INITIAL_MATRIX,&D)); 86 87 PetscCall(MatMatMult(A,D,MAT_INITIAL_MATRIX,1.0,&F)); 88 PetscCall(MatIsDiagonal(F)); 89 90 PetscCall(MatDestroy(&A)); 91 PetscCall(MatDestroy(&D)); 92 PetscCall(MatDestroy(&F)); 93 } 94 PetscCall(PetscFree(parts)); 95 PetscCall(PetscRandomDestroy(&rctx)); 96 97 PetscCall(PetscFinalize()); 98 return 0; 99 } 100 101 PetscErrorCode MatIsDiagonal(Mat A) 102 { 103 PetscInt ncols,i,j,rstart,rend; 104 const PetscInt *cols; 105 const PetscScalar *vals; 106 PetscBool founddiag; 107 108 PetscFunctionBeginUser; 109 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 110 for (i=rstart; i<rend; i++) { 111 founddiag = PETSC_FALSE; 112 PetscCall(MatGetRow(A,i,&ncols,&cols,&vals)); 113 for (j=0; j<ncols; j++) { 114 if (cols[j] == i) { 115 PetscCheck(PetscAbsScalar(vals[j] - 1) < PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row %" PetscInt_FMT " does not have 1 on the diagonal, it has %g",i,(double)PetscAbsScalar(vals[j])); 116 founddiag = PETSC_TRUE; 117 } else { 118 PetscCheck(PetscAbsScalar(vals[j]) < PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row %" PetscInt_FMT " has off-diagonal value %g at %" PetscInt_FMT "",i,(double)PetscAbsScalar(vals[j]),cols[j]); 119 } 120 } 121 PetscCall(MatRestoreRow(A,i,&ncols,&cols,&vals)); 122 PetscCheck(founddiag,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row %" PetscInt_FMT " does not have diagonal entrie",i); 123 } 124 PetscFunctionReturn(0); 125 } 126 127 /* 128 All processes receive all the block information 129 */ 130 PetscErrorCode BuildMatrix(const PetscInt *parts, PetscInt nblocks, const PetscInt *blocksizes,Mat *A) 131 { 132 PetscInt i,cnt = 0; 133 PetscMPIInt rank; 134 135 PetscFunctionBeginUser; 136 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 137 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD,parts[rank],parts[rank],PETSC_DETERMINE,PETSC_DETERMINE,0,NULL,0,NULL,A)); 138 PetscCall(MatSetOption(*A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 139 if (rank == 0) { 140 for (i=0; i<nblocks; i++) { 141 PetscCall(MatSetValue(*A,cnt,cnt+blocksizes[i]-1,1.0,INSERT_VALUES)); 142 PetscCall(MatSetValue(*A,cnt+blocksizes[i]-1,cnt,1.0,INSERT_VALUES)); 143 cnt += blocksizes[i]; 144 } 145 } 146 PetscCall(MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY)); 147 PetscCall(MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY)); 148 PetscCall(MatShift(*A,10)); 149 PetscFunctionReturn(0); 150 } 151 152 /*TEST 153 154 test: 155 156 test: 157 suffix: 2 158 nsize: 2 159 160 test: 161 suffix: 5 162 nsize: 5 163 164 TEST*/ 165