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