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