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