Actual source code: ex51.c
slepc-3.18.1 2022-11-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "Computes a partial GSVD of two matrices from IR Tools example.\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = number of grid subdivisions in x and y dimensions.\n\n";
15: #include <slepcsvd.h>
17: /* LookUp: returns an index i such that X(i) <= y < X(i+1), where X = linspace(0,1,N).
18: Only elements start..end-1 are considered */
19: PetscErrorCode LookUp(PetscInt N,PetscInt start,PetscInt end,PetscReal y,PetscInt *i)
20: {
21: PetscInt n=end-start,j=n/2;
22: PetscReal h=1.0/(N-1);
25: if (y<(start+j)*h) LookUp(N,start,start+j,y,i);
26: else if (y<(start+j+1)*h) *i = start+j;
27: else LookUp(N,start+j,end,y,i);
28: return 0;
29: }
31: /*
32: PermuteMatrices - Symmetric permutation of A using MatPartitioning, then permute
33: columns of B in the same way.
34: */
35: PetscErrorCode PermuteMatrices(Mat *A,Mat *B)
36: {
37: MatPartitioning part;
38: IS isn,is,id;
39: PetscInt *nlocal,Istart,Iend;
40: PetscMPIInt size,rank;
41: MPI_Comm comm;
42: Mat in=*A,out;
44: PetscObjectGetComm((PetscObject)in,&comm);
45: MPI_Comm_size(comm,&size);
46: MPI_Comm_rank(comm,&rank);
47: MatPartitioningCreate(comm,&part);
48: MatPartitioningSetAdjacency(part,in);
49: MatPartitioningSetFromOptions(part);
50: MatPartitioningApply(part,&is); /* get owner of each vertex */
51: ISPartitioningToNumbering(is,&isn); /* get new global number of each vertex */
52: PetscMalloc1(size,&nlocal);
53: ISPartitioningCount(is,size,nlocal); /* count vertices assigned to each process */
54: ISDestroy(&is);
56: /* get old global number of each new global number */
57: ISInvertPermutation(isn,nlocal[rank],&is);
58: PetscFree(nlocal);
59: ISDestroy(&isn);
60: MatPartitioningDestroy(&part);
62: /* symmetric permutation of A */
63: MatCreateSubMatrix(in,is,is,MAT_INITIAL_MATRIX,&out);
64: MatDestroy(A);
65: *A = out;
67: /* column permutation of B */
68: MatGetOwnershipRange(*B,&Istart,&Iend);
69: ISCreateStride(comm,Iend-Istart,Istart,1,&id);
70: ISSetIdentity(id);
71: MatCreateSubMatrix(*B,id,is,MAT_INITIAL_MATRIX,&out);
72: MatDestroy(B);
73: *B = out;
74: ISDestroy(&id);
75: ISDestroy(&is);
76: return 0;
77: }
79: int main(int argc,char **argv)
80: {
81: Mat A,B; /* operator matrices */
82: SVD svd; /* singular value problem solver context */
83: KSP ksp;
84: PetscInt n=32,N,i,i2,j,k,xidx,yidx,bl,Istart,Iend,col[3];
85: PetscScalar vals[] = { 1, -2, 1 },X,Y;
86: PetscBool flg,terse,permute=PETSC_FALSE;
87: PetscRandom rctx;
90: SlepcInitialize(&argc,&argv,(char*)0,help);
92: PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg);
93: N = n*n;
94: PetscPrintf(PETSC_COMM_WORLD,"\nGSVD of inverse interpolation problem, (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT "\n\n",N,2*N,N);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Build the matrices
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
101: PetscRandomSetInterval(rctx,0,1);
102: PetscRandomSetFromOptions(rctx);
104: MatCreate(PETSC_COMM_WORLD,&A);
105: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
106: MatSetFromOptions(A);
107: MatSetUp(A);
108: MatGetOwnershipRange(A,&Istart,&Iend);
110: /* make sure that the matrix is the same irrespective of the number of MPI processes */
111: PetscRandomSetSeed(rctx,0x12345678);
112: PetscRandomSeed(rctx);
113: for (k=0;k<Istart;k++) {
114: PetscRandomGetValue(rctx,&X);
115: PetscRandomGetValue(rctx,&Y);
116: }
118: for (k=0;k<Iend-Istart;k++) {
119: PetscRandomGetValue(rctx,&X);
120: LookUp(n,0,n,PetscRealPart(X),&xidx);
121: X = X*(n-1)-xidx; /* scale value to a 1-spaced grid */
122: PetscRandomGetValue(rctx,&Y);
123: LookUp(n,0,n,PetscRealPart(Y),&yidx);
124: Y = Y*(n-1)-yidx; /* scale value to a 1-spaced grid */
125: for (j=0;j<n;j++) {
126: for (i=0;i<n;i++) {
127: if (i<n-1 && j<n-1 && xidx==j && yidx==i) MatSetValue(A,Istart+k,i+j*n,1.0-X-Y+X*Y,ADD_VALUES);
128: if (i<n-1 && j>0 && xidx==j-1 && yidx==i) MatSetValue(A,Istart+k,i+j*n,X-X*Y,ADD_VALUES);
129: if (i>0 && j<n-1 && xidx==j && yidx==i-1) MatSetValue(A,Istart+k,i+j*n,Y-X*Y,ADD_VALUES);
130: if (i>0 && j>0 && xidx==j-1 && yidx==i-1) MatSetValue(A,Istart+k,i+j*n,X*Y,ADD_VALUES);
131: }
132: }
133: }
134: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
136: PetscRandomDestroy(&rctx);
138: MatCreate(PETSC_COMM_WORLD,&B);
139: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,2*N,N);
140: MatSetFromOptions(B);
141: MatSetUp(B);
143: for (i=Istart;i<Iend;i++) {
144: /* upper block: kron(speye(n),T1) where T1 is tridiagonal */
145: i2 = i+Istart;
146: if (i%n==0) MatSetValue(B,i2,i,1.0,INSERT_VALUES);
147: else if (i%n==n-1) {
148: MatSetValue(B,i2,i-1,-1.0,INSERT_VALUES);
149: MatSetValue(B,i2,i,1.0,INSERT_VALUES);
150: } else {
151: col[0]=i-1; col[1]=i; col[2]=i+1;
152: MatSetValues(B,1,&i2,3,col,vals,INSERT_VALUES);
153: }
154: /* lower block: kron(T2,speye(n)) where T2 is tridiagonal */
155: i2 = i+Iend;
156: bl = i/n; /* index of block */
157: j = i-bl*n; /* index within block */
158: if (bl==0 || bl==n-1) MatSetValue(B,i2,i,1.0,INSERT_VALUES);
159: else {
160: col[0]=i-n; col[1]=i; col[2]=i+n;
161: MatSetValues(B,1,&i2,3,col,vals,INSERT_VALUES);
162: }
163: }
164: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
165: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
167: PetscOptionsGetBool(NULL,NULL,"-permute",&permute,NULL);
168: if (permute) {
169: PetscPrintf(PETSC_COMM_WORLD," Permuting to optimize parallel matvec\n");
170: PermuteMatrices(&A,&B);
171: }
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Create the singular value solver and set various options
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: SVDCreate(PETSC_COMM_WORLD,&svd);
178: SVDSetOperators(svd,A,B);
179: SVDSetProblemType(svd,SVD_GENERALIZED);
181: SVDSetType(svd,SVDTRLANCZOS);
182: SVDSetDimensions(svd,6,PETSC_DEFAULT,PETSC_DEFAULT);
183: SVDTRLanczosSetExplicitMatrix(svd,PETSC_TRUE);
184: SVDTRLanczosSetScale(svd,-10);
185: SVDTRLanczosGetKSP(svd,&ksp);
186: KSPSetTolerances(ksp,1e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
188: SVDSetFromOptions(svd);
190: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191: Solve the problem and print solution
192: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194: SVDSolve(svd);
196: /* show detailed info unless -terse option is given by user */
197: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
198: if (terse) SVDErrorView(svd,SVD_ERROR_NORM,NULL);
199: else {
200: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
201: SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD);
202: SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD);
203: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
204: }
205: SVDDestroy(&svd);
206: MatDestroy(&A);
207: MatDestroy(&B);
208: SlepcFinalize();
209: return 0;
210: }
212: /*TEST
214: testset:
215: args: -n 16 -terse
216: requires: double
217: output_file: output/ex51_1.out
218: test:
219: args: -svd_trlanczos_gbidiag {{upper lower}} -svd_trlanczos_oneside
220: suffix: 1
221: test:
222: suffix: 2
223: nsize: 2
224: args: -permute
225: filter: grep -v Permuting
227: TEST*/