Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Tests multiple calls to SVDSolve with equal matrix size (GSVD).\n\n"
12 : "The command line options are:\n"
13 : " -m <m>, where <m> = number of rows of A.\n"
14 : " -n <n>, where <n> = number of columns of A.\n"
15 : " -p <p>, where <p> = number of rows of B.\n\n";
16 :
17 : #include <slepcsvd.h>
18 :
19 : /*
20 : This example solves two GSVD problems for the bidiagonal matrices
21 :
22 : | 1 2 | | 1 |
23 : | 1 2 | | 2 1 |
24 : | 1 2 | | 2 1 |
25 : A1 = | . . | A2 = | . . |
26 : | . . | | . . |
27 : | 1 2 | | 2 1 |
28 : | 1 2 | | 2 1 |
29 :
30 : with B = tril(ones(p,n))
31 : */
32 :
33 11 : int main(int argc,char **argv)
34 : {
35 11 : Mat A1,A2,B;
36 11 : SVD svd;
37 11 : PetscInt m=15,n=20,p=21,Istart,Iend,i,j,d,col[2];
38 11 : PetscScalar valsa[] = { 1, 2 }, valsb[] = { 2, 1 };
39 :
40 11 : PetscFunctionBeginUser;
41 11 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
42 11 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
43 11 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
44 11 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL));
45 11 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized singular value decomposition, (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT "\n\n",m,p,n));
46 :
47 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48 : Generate the matrices
49 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50 :
51 11 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A1));
52 11 : PetscCall(MatSetSizes(A1,PETSC_DECIDE,PETSC_DECIDE,m,n));
53 11 : PetscCall(MatSetFromOptions(A1));
54 11 : PetscCall(MatGetOwnershipRange(A1,&Istart,&Iend));
55 146 : for (i=Istart;i<Iend;i++) {
56 135 : col[0]=i; col[1]=i+1;
57 135 : if (i<n-1) PetscCall(MatSetValues(A1,1,&i,2,col,valsa,INSERT_VALUES));
58 135 : else if (i==n-1) PetscCall(MatSetValue(A1,i,col[0],valsa[0],INSERT_VALUES));
59 : }
60 11 : PetscCall(MatAssemblyBegin(A1,MAT_FINAL_ASSEMBLY));
61 11 : PetscCall(MatAssemblyEnd(A1,MAT_FINAL_ASSEMBLY));
62 :
63 11 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A2));
64 11 : PetscCall(MatSetSizes(A2,PETSC_DECIDE,PETSC_DECIDE,m,n));
65 11 : PetscCall(MatSetFromOptions(A2));
66 11 : PetscCall(MatGetOwnershipRange(A2,&Istart,&Iend));
67 146 : for (i=Istart;i<Iend;i++) {
68 135 : col[0]=i-1; col[1]=i;
69 135 : if (i==0) PetscCall(MatSetValue(A2,i,col[1],valsb[1],INSERT_VALUES));
70 126 : else if (i<n) PetscCall(MatSetValues(A2,1,&i,2,col,valsb,INSERT_VALUES));
71 135 : else if (i==n) PetscCall(MatSetValue(A2,i,col[0],valsb[0],INSERT_VALUES));
72 : }
73 11 : PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY));
74 11 : PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY));
75 :
76 11 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
77 11 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
78 11 : PetscCall(MatSetFromOptions(B));
79 11 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
80 11 : d = PetscMax(0,n-p);
81 200 : for (i=Istart;i<Iend;i++) {
82 1179 : for (j=PetscMax(0,i-5);j<=PetscMin(i,n-1);j++) PetscCall(MatSetValue(B,i,j+d,1.0,INSERT_VALUES));
83 : }
84 11 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
85 11 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
86 :
87 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88 : Create the singular value solver, set options and solve
89 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90 :
91 11 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
92 11 : PetscCall(SVDSetOperators(svd,A1,B));
93 11 : PetscCall(SVDSetFromOptions(svd));
94 11 : PetscCall(SVDSolve(svd));
95 11 : PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
96 :
97 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98 : Solve second problem
99 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 :
101 11 : PetscCall(SVDSetOperators(svd,A2,B));
102 11 : PetscCall(SVDSolve(svd));
103 11 : PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
104 :
105 : /* Free work space */
106 11 : PetscCall(SVDDestroy(&svd));
107 11 : PetscCall(MatDestroy(&A1));
108 11 : PetscCall(MatDestroy(&A2));
109 11 : PetscCall(MatDestroy(&B));
110 11 : PetscCall(SlepcFinalize());
111 : return 0;
112 : }
113 :
114 : /*TEST
115 :
116 : testset:
117 : args: -svd_nsv 3
118 : requires: !single
119 : output_file: output/test16_1.out
120 : test:
121 : suffix: 1_lapack
122 : args: -svd_type lapack
123 : test:
124 : suffix: 1_cross
125 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
126 : test:
127 : suffix: 1_cyclic
128 : args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
129 : test:
130 : suffix: 1_trlanczos
131 : args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single lower}} -svd_trlanczos_ksp_rtol 1e-10
132 : requires: double
133 : test:
134 : suffix: 1_trlanczos_par
135 : nsize: 2
136 : args: -svd_type trlanczos -ds_parallel {{redundant synchronized}}
137 :
138 : testset:
139 : args: -svd_nsv 3 -mat_type aijcusparse
140 : requires: cuda !single
141 : output_file: output/test16_1.out
142 : test:
143 : suffix: 2_cross
144 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
145 : test:
146 : suffix: 2_cyclic
147 : args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
148 : test:
149 : suffix: 2_trlanczos
150 : args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single lower}} -svd_trlanczos_ksp_rtol 1e-10
151 : requires: double
152 :
153 : TEST*/
|
