1 static char help[] = "Test ViennaCL Matrix Conversions"; 2 3 #include <petscmat.h> 4 5 int main(int argc,char **args) 6 { 7 PetscMPIInt rank,size; 8 9 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 10 11 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 12 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 13 14 /* Construct a sequential ViennaCL matrix and vector */ 15 if (rank == 0) { 16 Mat A_vcl; 17 Vec v_vcl,r_vcl; 18 PetscInt n = 17, m = 31,nz = 5,i,cnt,j; 19 const PetscScalar val = 1.0; 20 21 CHKERRQ(MatCreateSeqAIJViennaCL(PETSC_COMM_SELF,m,n,nz,NULL,&A_vcl)); 22 23 /* Add nz arbitrary entries per row in arbitrary columns */ 24 for (i=0;i<m;++i) { 25 for (cnt = 0; cnt<nz; ++cnt) { 26 j = (19 * cnt + (7*i + 3)) % n; 27 CHKERRQ(MatSetValue(A_vcl,i,j,(PetscScalar)(0.3 * i + j),INSERT_VALUES)); 28 } 29 } 30 CHKERRQ(MatAssemblyBegin(A_vcl,MAT_FINAL_ASSEMBLY)); 31 CHKERRQ(MatAssemblyEnd(A_vcl,MAT_FINAL_ASSEMBLY)); 32 33 CHKERRQ(VecCreateSeqViennaCL(PETSC_COMM_SELF,n,&v_vcl)); 34 CHKERRQ(VecCreateSeqViennaCL(PETSC_COMM_SELF,m,&r_vcl)); 35 CHKERRQ(VecSet(v_vcl,val)); 36 37 CHKERRQ(MatMult(A_vcl,v_vcl,r_vcl)); 38 39 CHKERRQ(VecDestroy(&v_vcl)); 40 CHKERRQ(VecDestroy(&r_vcl)); 41 CHKERRQ(MatDestroy(&A_vcl)); 42 } 43 44 /* Create a sequential AIJ matrix on rank 0 convert it to a new ViennaCL matrix, and apply it to a ViennaCL vector */ 45 if (rank == 0) { 46 Mat A,A_vcl; 47 Vec v,r,v_vcl,r_vcl,d_vcl; 48 PetscInt n=17,m=31,nz=5,i,cnt,j; 49 const PetscScalar val = 1.0; 50 PetscReal dnorm; 51 const PetscReal tol = 1e-5; 52 53 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,m,n,nz,NULL,&A)); 54 55 /* Add nz arbitrary entries per row in arbitrary columns */ 56 for (i=0;i<m;++i) { 57 for (cnt = 0; cnt<nz; ++cnt) { 58 j = (19 * cnt + (7*i + 3)) % n; 59 CHKERRQ(MatSetValue(A,i,j,(PetscScalar) (0.3 * i + j),INSERT_VALUES)); 60 } 61 } 62 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 63 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 64 65 CHKERRQ(PetscObjectSetName((PetscObject)A,"Sequential CPU Matrix")); 66 67 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&v)); 68 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,m,&r)); 69 CHKERRQ(PetscObjectSetName((PetscObject)r,"CPU result vector")); 70 CHKERRQ(VecSet(v,val)); 71 CHKERRQ(MatMult(A,v,r)); 72 73 CHKERRQ(MatConvert(A,MATSEQAIJVIENNACL,MAT_INITIAL_MATRIX,&A_vcl)); 74 CHKERRQ(PetscObjectSetName((PetscObject)A_vcl,"New ViennaCL Matrix")); 75 76 CHKERRQ(VecCreateSeqViennaCL(PETSC_COMM_SELF,n,&v_vcl)); 77 CHKERRQ(VecCreateSeqViennaCL(PETSC_COMM_SELF,m,&r_vcl)); 78 CHKERRQ(PetscObjectSetName((PetscObject)r_vcl,"ViennaCL result vector")); 79 CHKERRQ(VecSet(v_vcl,val)); 80 CHKERRQ(MatMult(A_vcl,v_vcl,r_vcl)); 81 82 CHKERRQ(VecDuplicate(r_vcl,&d_vcl)); 83 CHKERRQ(VecCopy(r_vcl,d_vcl)); 84 CHKERRQ(VecAXPY(d_vcl,-1.0,r_vcl)); 85 CHKERRQ(VecNorm(d_vcl,NORM_INFINITY,&dnorm)); 86 PetscCheckFalse(dnorm > tol,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"Sequential CPU and MPI ViennaCL vector results incompatible with inf norm greater than tolerance of %g",tol); 87 88 CHKERRQ(VecDestroy(&v)); 89 CHKERRQ(VecDestroy(&r)); 90 CHKERRQ(VecDestroy(&v_vcl)); 91 CHKERRQ(VecDestroy(&r_vcl)); 92 CHKERRQ(VecDestroy(&d_vcl)); 93 CHKERRQ(MatDestroy(&A)); 94 CHKERRQ(MatDestroy(&A_vcl)); 95 } 96 97 /* Create a sequential AIJ matrix on rank 0 convert it inplace to a new ViennaCL matrix, and apply it to a ViennaCL vector */ 98 if (rank == 0) { 99 Mat A; 100 Vec v,r,v_vcl,r_vcl,d_vcl; 101 PetscInt n=17,m=31,nz=5,i,cnt,j; 102 const PetscScalar val = 1.0; 103 PetscReal dnorm; 104 const PetscReal tol = 1e-5; 105 106 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,m,n,nz,NULL,&A)); 107 108 /* Add nz arbitrary entries per row in arbitrary columns */ 109 for (i=0;i<m;++i) { 110 for (cnt = 0; cnt<nz; ++cnt) { 111 j = (19 * cnt + (7*i + 3)) % n; 112 CHKERRQ(MatSetValue(A,i,j,(PetscScalar)(0.3 * i + j),INSERT_VALUES)); 113 } 114 } 115 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 116 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 117 118 CHKERRQ(PetscObjectSetName((PetscObject)A,"Sequential CPU Matrix")); 119 120 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&v)); 121 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,m,&r)); 122 CHKERRQ(PetscObjectSetName((PetscObject)r,"CPU result vector")); 123 CHKERRQ(VecSet(v,val)); 124 CHKERRQ(MatMult(A,v,r)); 125 126 CHKERRQ(MatConvert(A,MATSEQAIJVIENNACL,MAT_INPLACE_MATRIX,&A)); 127 CHKERRQ(PetscObjectSetName((PetscObject)A,"Converted ViennaCL Matrix")); 128 129 CHKERRQ(VecCreateSeqViennaCL(PETSC_COMM_SELF,n,&v_vcl)); 130 CHKERRQ(VecCreateSeqViennaCL(PETSC_COMM_SELF,m,&r_vcl)); 131 CHKERRQ(PetscObjectSetName((PetscObject)r_vcl,"ViennaCL result vector")); 132 CHKERRQ(VecSet(v_vcl,val)); 133 CHKERRQ(MatMult(A,v_vcl,r_vcl)); 134 135 CHKERRQ(VecDuplicate(r_vcl,&d_vcl)); 136 CHKERRQ(VecCopy(r_vcl,d_vcl)); 137 CHKERRQ(VecAXPY(d_vcl,-1.0,r_vcl)); 138 CHKERRQ(VecNorm(d_vcl,NORM_INFINITY,&dnorm)); 139 PetscCheckFalse(dnorm > tol,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MPI CPU and MPI ViennaCL Vector results incompatible with inf norm greater than tolerance of %g",tol); 140 141 CHKERRQ(VecDestroy(&v)); 142 CHKERRQ(VecDestroy(&r)); 143 CHKERRQ(VecDestroy(&v_vcl)); 144 CHKERRQ(VecDestroy(&r_vcl)); 145 CHKERRQ(VecDestroy(&d_vcl)); 146 CHKERRQ(MatDestroy(&A)); 147 } 148 149 /* Create a parallel AIJ matrix, convert it to a new ViennaCL matrix, and apply it to a ViennaCL vector */ 150 if (size > 1) { 151 Mat A,A_vcl; 152 Vec v,r,v_vcl,r_vcl,d_vcl; 153 PetscInt N=17,M=31,nz=5,i,cnt,j,rlow,rhigh; 154 const PetscScalar val = 1.0; 155 PetscReal dnorm; 156 const PetscReal tol=1e-5; 157 158 CHKERRQ(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,M,N,nz,NULL,nz,NULL,&A)); 159 160 /* Add nz arbitrary entries per row in arbitrary columns */ 161 CHKERRQ(MatGetOwnershipRange(A,&rlow,&rhigh)); 162 for (i=rlow;i<rhigh;++i) { 163 for (cnt = 0; cnt<nz; ++cnt) { 164 j = (19 * cnt + (7*i + 3)) % N; 165 CHKERRQ(MatSetValue(A,i,j,(PetscScalar)(0.3 * i + j),INSERT_VALUES)); 166 } 167 } 168 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 169 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 170 171 CHKERRQ(PetscObjectSetName((PetscObject)A,"MPI CPU Matrix")); 172 173 CHKERRQ(MatCreateVecs(A,&v,&r)); 174 CHKERRQ(PetscObjectSetName((PetscObject)r,"MPI CPU result vector")); 175 CHKERRQ(VecSet(v,val)); 176 CHKERRQ(MatMult(A,v,r)); 177 178 CHKERRQ(MatConvert(A,MATMPIAIJVIENNACL,MAT_INITIAL_MATRIX,&A_vcl)); 179 CHKERRQ(PetscObjectSetName((PetscObject)A_vcl,"New MPI ViennaCL Matrix")); 180 181 CHKERRQ(MatCreateVecs(A_vcl,&v_vcl,&r_vcl)); 182 CHKERRQ(PetscObjectSetName((PetscObject)r_vcl,"ViennaCL result vector")); 183 CHKERRQ(VecSet(v_vcl,val)); 184 CHKERRQ(MatMult(A_vcl,v_vcl,r_vcl)); 185 186 CHKERRQ(VecDuplicate(r_vcl,&d_vcl)); 187 CHKERRQ(VecCopy(r_vcl,d_vcl)); 188 CHKERRQ(VecAXPY(d_vcl,-1.0,r_vcl)); 189 CHKERRQ(VecNorm(d_vcl,NORM_INFINITY,&dnorm)); 190 PetscCheckFalse(dnorm > tol,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MPI CPU and MPI ViennaCL Vector results incompatible with inf norm greater than tolerance of %g",tol); 191 192 CHKERRQ(VecDestroy(&v)); 193 CHKERRQ(VecDestroy(&r)); 194 CHKERRQ(VecDestroy(&v_vcl)); 195 CHKERRQ(VecDestroy(&r_vcl)); 196 CHKERRQ(VecDestroy(&d_vcl)); 197 CHKERRQ(MatDestroy(&A)); 198 CHKERRQ(MatDestroy(&A_vcl)); 199 } 200 201 /* Create a parallel AIJ matrix, convert it in-place to a ViennaCL matrix, and apply it to a ViennaCL vector */ 202 if (size > 1) { 203 Mat A; 204 Vec v,r,v_vcl,r_vcl,d_vcl; 205 PetscInt N=17,M=31,nz=5,i,cnt,j,rlow,rhigh; 206 const PetscScalar val = 1.0; 207 PetscReal dnorm; 208 const PetscReal tol=1e-5; 209 210 CHKERRQ(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,M,N,nz,NULL,nz,NULL,&A)); 211 212 /* Add nz arbitrary entries per row in arbitrary columns */ 213 CHKERRQ(MatGetOwnershipRange(A,&rlow,&rhigh)); 214 for (i=rlow;i<rhigh;++i) { 215 for (cnt = 0; cnt<nz; ++cnt) { 216 j = (19 * cnt + (7*i + 3)) % N; 217 CHKERRQ(MatSetValue(A,i,j,(PetscScalar)(0.3 * i + j),INSERT_VALUES)); 218 } 219 } 220 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 221 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 222 223 CHKERRQ(PetscObjectSetName((PetscObject)A,"MPI CPU Matrix")); 224 225 CHKERRQ(MatCreateVecs(A,&v,&r)); 226 CHKERRQ(PetscObjectSetName((PetscObject)r,"MPI CPU result vector")); 227 CHKERRQ(VecSet(v,val)); 228 CHKERRQ(MatMult(A,v,r)); 229 230 CHKERRQ(MatConvert(A,MATMPIAIJVIENNACL,MAT_INPLACE_MATRIX,&A)); 231 CHKERRQ(PetscObjectSetName((PetscObject)A,"Converted MPI ViennaCL Matrix")); 232 233 CHKERRQ(MatCreateVecs(A,&v_vcl,&r_vcl)); 234 CHKERRQ(PetscObjectSetName((PetscObject)r_vcl,"ViennaCL result vector")); 235 CHKERRQ(VecSet(v_vcl,val)); 236 CHKERRQ(MatMult(A,v_vcl,r_vcl)); 237 238 CHKERRQ(VecDuplicate(r_vcl,&d_vcl)); 239 CHKERRQ(VecCopy(r_vcl,d_vcl)); 240 CHKERRQ(VecAXPY(d_vcl,-1.0,r_vcl)); 241 CHKERRQ(VecNorm(d_vcl,NORM_INFINITY,&dnorm)); 242 PetscCheckFalse(dnorm > tol,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MPI CPU and MPI ViennaCL Vector results incompatible with inf norm greater than tolerance of %g",tol); 243 244 CHKERRQ(VecDestroy(&v)); 245 CHKERRQ(VecDestroy(&r)); 246 CHKERRQ(VecDestroy(&v_vcl)); 247 CHKERRQ(VecDestroy(&r_vcl)); 248 CHKERRQ(VecDestroy(&d_vcl)); 249 CHKERRQ(MatDestroy(&A)); 250 } 251 252 CHKERRQ(PetscFinalize()); 253 return 0; 254 255 } 256 257 /*TEST 258 259 build: 260 requires: viennacl 261 262 test: 263 output_file: output/ex204.out 264 265 test: 266 suffix: 2 267 nsize: 2 268 output_file: output/ex204.out 269 270 TEST*/ 271