1 static char help[] = "Test MatMult() for Hermitian matrix.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc,char **args) 6 { 7 Mat A,As; 8 PetscBool flg; 9 PetscMPIInt size; 10 PetscInt i,j; 11 PetscScalar v,sigma2; 12 PetscReal h2,sigma1=100.0; 13 PetscInt dim,Ii,J,n = 3,rstart,rend; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 17 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 18 PetscCall(PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL)); 19 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 20 dim = n*n; 21 22 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 23 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim)); 24 PetscCall(MatSetType(A,MATAIJ)); 25 PetscCall(MatSetFromOptions(A)); 26 PetscCall(MatSetUp(A)); 27 28 sigma2 = 10.0*PETSC_i; 29 h2 = 1.0/((n+1)*(n+1)); 30 31 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 32 for (Ii=rstart; Ii<rend; Ii++) { 33 v = -1.0; i = Ii/n; j = Ii - i*n; 34 if (i>0) { 35 J = Ii-n; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 36 } 37 if (i<n-1) { 38 J = Ii+n; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 39 } 40 if (j>0) { 41 J = Ii-1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 42 } 43 if (j<n-1) { 44 J = Ii+1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 45 } 46 v = 4.0 - sigma1*h2; 47 PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES)); 48 } 49 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 50 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 51 52 /* Check whether A is symmetric */ 53 PetscCall(PetscOptionsHasName(NULL,NULL, "-check_symmetric", &flg)); 54 if (flg) { 55 PetscCall(MatIsSymmetric(A,0.0,&flg)); 56 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_USER,"A is not symmetric"); 57 } 58 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 59 60 /* make A complex Hermitian */ 61 Ii = 0; J = dim-1; 62 if (Ii >= rstart && Ii < rend) { 63 v = sigma2*h2; /* RealPart(v) = 0.0 */ 64 PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 65 v = -sigma2*h2; 66 PetscCall(MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES)); 67 } 68 69 Ii = dim-2; J = dim-1; 70 if (Ii >= rstart && Ii < rend) { 71 v = sigma2*h2; /* RealPart(v) = 0.0 */ 72 PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 73 v = -sigma2*h2; 74 PetscCall(MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES)); 75 } 76 77 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 78 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 79 PetscCall(MatViewFromOptions(A,NULL,"-disp_mat")); 80 81 /* Check whether A is Hermitian, then set A->hermitian flag */ 82 PetscCall(PetscOptionsHasName(NULL,NULL, "-check_Hermitian", &flg)); 83 if (flg && size == 1) { 84 PetscCall(MatIsHermitian(A,0.0,&flg)); 85 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_USER,"A is not Hermitian"); 86 } 87 PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE)); 88 89 #if defined(PETSC_HAVE_SUPERLU_DIST) 90 /* Test Cholesky factorization */ 91 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_choleskyfactor", &flg)); 92 if (flg) { 93 Mat F; 94 IS perm,iperm; 95 MatFactorInfo info; 96 PetscInt nneg,nzero,npos; 97 98 PetscCall(MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_CHOLESKY,&F)); 99 PetscCall(MatGetOrdering(A,MATORDERINGND,&perm,&iperm)); 100 PetscCall(MatCholeskyFactorSymbolic(F,A,perm,&info)); 101 PetscCall(MatCholeskyFactorNumeric(F,A,&info)); 102 103 PetscCall(MatGetInertia(F,&nneg,&nzero,&npos)); 104 PetscCall(PetscPrintf(PETSC_COMM_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos)); 105 PetscCall(MatDestroy(&F)); 106 PetscCall(ISDestroy(&perm)); 107 PetscCall(ISDestroy(&iperm)); 108 } 109 #endif 110 111 /* Create a Hermitian matrix As in sbaij format */ 112 PetscCall(MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As)); 113 PetscCall(MatViewFromOptions(As,NULL,"-disp_mat")); 114 115 /* Test MatMult */ 116 PetscCall(MatMultEqual(A,As,10,&flg)); 117 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatMult not equal"); 118 PetscCall(MatMultAddEqual(A,As,10,&flg)); 119 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatMultAdd not equal"); 120 121 /* Free spaces */ 122 PetscCall(MatDestroy(&A)); 123 PetscCall(MatDestroy(&As)); 124 PetscCall(PetscFinalize()); 125 return 0; 126 } 127 128 /*TEST 129 130 build: 131 requires: complex 132 133 test: 134 args: -n 1000 135 output_file: output/ex127.out 136 137 test: 138 suffix: 2 139 nsize: 3 140 args: -n 1000 141 output_file: output/ex127.out 142 143 test: 144 suffix: superlu_dist 145 nsize: 3 146 requires: superlu_dist 147 args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM 148 output_file: output/ex127_superlu_dist.out 149 TEST*/ 150