Actual source code: ex53.c
slepc-main 2024-12-17
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[] = "Partial hyperbolic singular value decomposition (HSVD) of matrix generated by plane rotations.\n"
12: "The command line options are:\n"
13: " -p <p>, number of rows in upper part of matrix A.\n"
14: " -q <q>, number of rows in lower part of matrix A.\n"
15: " -n <n>, number of columns of matrix A.\n"
16: " -k <k>, condition number.\n"
17: " -d <d>, density.\n"
18: " -prog, show progress of matrix generation\n";
20: #include <slepcsvd.h>
21: #include <slepcblaslapack.h>
23: /* The problem matrix is the (p+q)-by-n matrix
25: A = | Q1 * D1 * U |
26: | 1/2*Q2 * D2 * U |
28: with p>=n, where Q1 (p-by-p), Q2 (q-by-q) and U (n-by-n) are orthonormal
29: matrices, D1 is a diagonal p-by-n matrix with elements
30: D1(i,i) = sqrt(5/4)*d(i), if i<=q
31: = d(i), if i>q
32: D2 is a diagonal q-by-n matrix with elements D2(i,i)= d(i),
33: and d is an n-element vector with elements distributed exponentially
34: from 1/k to 1.
36: The signature matrix is Omega = diag(I_p,-I_q)
38: The hyperbolic singular values of A are equal to the elements of vector d.
40: A is generated by random plane rotations applied to the diagonal matrices
41: D1 and D2
43: This test case is based on the paper:
45: [1] A. Bojanczyk, N.J. Higham, and H. Patel, "Solving the Indefinite Least
46: Squares Problem by Hyperbolic QR Factorization", SIAM J. Matrix Anal.
47: Appl. 24(4):914-931 2003
48: */
50: /* Timing */
51: #ifdef TIMING
52: static PetscReal tt_rotr,tt_rotc,tt_assmr=0.0,tt_assmc=0.0,tt_getr=0.0,tt_getc=0.0,tt_setvalr=0.0,tt_setvalc=0.0;
53: #define time_now(t) t=MPI_Wtime();
54: #define time_diff(tacum,t0,t,t1) {t=MPI_Wtime(); tacum += t-t0; t1=t;}
55: #else
56: #define time_diff(tacum,t0,t,t1)
57: #define time_now(t)
58: #endif
60: struct _p_XMat {
61: MPI_Comm comm;
62: PetscInt M,N; /* global number of rows and columns */
63: PetscInt Istart,Iend; /* Index of 1st row and last local rows */
64: PetscInt *nzr; /* number of non-zeros in each row */
65: PetscInt **rcol; /* row columns */
66: PetscInt *ranges; /* proc row ranges */
67: PetscScalar **rval; /* row values */
68: };
70: typedef struct _p_XMat *XMat;
72: /* Iterator over non-zero rows of 2 columns */
73: struct _p_ColsNzIter {
74: XMat A;
75: PetscInt j1,j2; /* Columns considered */
76: PetscInt iloc; /* current local row */
77: };
79: typedef struct _p_ColsNzIter *ColsNzIter;
81: static PetscErrorCode XMatCreate(MPI_Comm comm,XMat *A,PetscInt Istart,PetscInt Iend,PetscInt M,PetscInt N)
82: {
83: PetscInt i;
84: PetscMPIInt np;
86: PetscFunctionBeginUser;
87: PetscCall(PetscNew(A));
88: (*A)->M = M;
89: (*A)->N = N;
90: (*A)->Istart = Istart;
91: (*A)->Iend = Iend;
92: (*A)->comm = comm;
93: PetscCallMPI(MPI_Comm_size(comm,&np));
94: PetscCall(PetscMalloc4(Iend-Istart,&(*A)->nzr,Iend-Istart,&(*A)->rcol,Iend-Istart,&(*A)->rval,np+1,&(*A)->ranges));
95: for (i=0; i<Iend-Istart; i++) {
96: (*A)->nzr[i] = 0;
97: (*A)->rcol[i] = NULL;
98: (*A)->rval[i] = NULL;
99: }
100: PetscCallMPI(MPI_Allgather(&Istart,1,MPIU_INT,(*A)->ranges,1,MPIU_INT,comm));
101: (*A)->ranges[np]=M;
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode XMatDestroy(XMat *A)
106: {
107: PetscInt m,i;
109: PetscFunctionBeginUser;
110: m = (*A)->Iend-(*A)->Istart;
111: for (i=0; i<m; i++) {
112: if ((*A)->nzr[i]>0) {
113: PetscCall(PetscFree2((*A)->rcol[i],(*A)->rval[i]));
114: }
115: }
116: PetscCall(PetscFree4((*A)->nzr,(*A)->rcol,(*A)->rval,(*A)->ranges));
117: PetscCall(PetscFree(*A));
118: *A = NULL;
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode XMatGetSize(XMat A,PetscInt *M,PetscInt *N)
123: {
124: PetscFunctionBeginUser;
125: *M = A->M;
126: *N = A->N;
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode XMatGetOwnershipRange(XMat A,PetscInt *Istart,PetscInt *Iend)
131: {
132: PetscFunctionBeginUser;
133: *Istart = A->Istart;
134: *Iend = A->Iend;
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode XMatGetOwnershipRanges(XMat A,PetscInt **ranges)
139: {
140: PetscFunctionBeginUser;
141: *ranges = A->ranges;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode XMatGetRow(XMat A,PetscInt i,PetscInt *nc,PetscInt **vj,PetscScalar **va)
146: {
147: PetscInt iloc=i-A->Istart;
149: PetscFunctionBeginUser;
150: *nc = A->nzr[iloc];
151: *vj = A->rcol[iloc];
152: *va = A->rval[iloc];
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode XMatRestoreRow(XMat A,PetscInt i,PetscInt *nc,PetscInt **vj,PetscScalar **va)
157: {
158: PetscFunctionBeginUser;
159: *nc = 0;
160: *vj = NULL;
161: *va = NULL;
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode XMatChangeRow(XMat A,PetscInt i,PetscInt nc,PetscInt vj[],PetscScalar va[])
166: {
167: PetscInt iloc = i-A->Istart;
169: PetscFunctionBeginUser;
170: if (A->nzr[iloc]>0) {
171: PetscCall(PetscFree2(A->rcol[iloc],A->rval[iloc]));
172: }
173: A->nzr[iloc] = nc;
174: A->rcol[iloc] = vj;
175: A->rval[iloc] = va;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /* Given a row with non-zero elements in columns indicated by rcol, get the
180: position where the non-zero element of column j is, or where it should be inserted.
181: Returns true if there is a non-zero element for the column, false otherwise */
182: static PetscBool getcolpos(const PetscInt rcol[],PetscInt nz,PetscInt j,PetscInt *k)
183: {
184: PetscInt ka,kb;
185: if (nz==0) {
186: *k = 0;
187: return PETSC_FALSE;
188: }
189: if (j<=rcol[0]) {
190: *k = 0;
191: return (rcol[0]==j)? PETSC_TRUE: PETSC_FALSE;
192: }
193: if (rcol[nz-1]<=j) {
194: *k = nz-(rcol[nz-1]==j);
195: return (rcol[nz-1]==j)? PETSC_TRUE: PETSC_FALSE;
196: }
197: /* Bisection: rcol[ka] < j < rcol[kb] */
198: ka = 0; kb = nz-1;
199: while (ka+1<kb) {
200: *k = (ka+kb)/2;
201: if (rcol[*k]<=j) {
202: if (rcol[*k]==j) return PETSC_TRUE;
203: else ka = *k;
204: } else kb = *k;
205: }
206: *k = kb;
207: return PETSC_FALSE;
208: }
210: /* Only local elements can be set */
211: static PetscErrorCode XMatSetValue(XMat A,PetscInt i,PetscInt j,PetscScalar x,InsertMode dummy)
212: {
213: PetscInt nz,iloc,k,kx,*col1,*col2;
214: PetscScalar *val1,*val2;
216: PetscFunctionBeginUser;
217: PetscCheck(i>=A->Istart && i<A->Iend && j>=0 && j<A->N,A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Incorrect index of row/column");
218: iloc = i-A->Istart;
219: nz = A->nzr[iloc];
220: col1 = A->rcol[iloc];
221: val1 = A->rval[iloc];
222: if (!getcolpos(col1,nz,j,&k)) {
223: A->nzr[iloc]++;
224: /* Replace row */
225: PetscCall(PetscMalloc2(nz+1,&col2,nz+1,&val2));
226: A->rcol[iloc] = col2;
227: A->rval[iloc] = val2;
228: for (kx=0; kx<k; kx++) {
229: col2[kx] = col1[kx];
230: val2[kx] = val1[kx];
231: }
232: for (; kx<nz; kx++) {
233: col2[kx+1] = col1[kx];
234: val2[kx+1] = val1[kx];
235: }
236: if (nz>0) {
237: PetscCall(PetscFree2(col1,val1));
238: }
239: col1 = col2;
240: val1 = val2;
241: }
242: col1[k] = j;
243: val1[k] = x;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /* Convert to PETSc Mat */
248: static PetscErrorCode XMatConvert(XMat XA,Mat A)
249: {
250: PetscInt i,iloc;
252: PetscFunctionBeginUser;
253: for (i=XA->Istart; i<XA->Iend; i++) {
254: iloc = i-XA->Istart;
255: PetscCall(MatSetValues(A,1,&i,XA->nzr[iloc],XA->rcol[iloc],XA->rval[iloc],INSERT_VALUES));
256: }
257: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
258: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /* gets the number of nonzeros in the matrix on this MPI rank */
263: static PetscErrorCode XMatGetNumberNonzeros(XMat A,PetscInt *nz)
264: {
265: PetscInt i;
267: PetscFunctionBeginUser;
268: *nz = 0;
269: for (i=0; i<A->Iend-A->Istart; i++) *nz += A->nzr[i];
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode XMatCreateColsNzIter(XMat A,PetscInt j1,PetscInt j2,ColsNzIter *iter)
274: {
275: PetscFunctionBeginUser;
276: PetscCall(PetscNew(iter));
277: (*iter)->A = A;
278: (*iter)->j1 = j1;
279: (*iter)->j2 = j2;
280: (*iter)->iloc = 0;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode ColitDestroy(ColsNzIter *iter)
285: {
286: PetscFunctionBeginUser;
287: PetscCall(PetscFree(*iter));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode ColitNextRow(ColsNzIter iter,PetscInt *i,PetscScalar **paij1,PetscScalar **paij2)
292: {
293: PetscInt iloc = iter->iloc,*rcol,nz,k,k1,k2,kx,jx,*rcol2;
294: PetscScalar *rval,*rval2;
295: PetscBool found1, found2;
296: XMat A=iter->A;
298: PetscFunctionBeginUser;
299: *i = -1;
300: *paij1 = *paij2 = NULL;
301: if (iloc>=0) {
302: rcol = A->rcol[iloc];
303: nz = A->nzr[iloc];
304: while (PETSC_TRUE) {
305: found1 = getcolpos(rcol,nz,iter->j1,&k1);
306: found2 = getcolpos(PetscSafePointerPlusOffset(rcol,k1+found1),nz-k1-found1,iter->j2,&k2);
307: if (found1 || found2) break;
308: iloc++;
309: if (iloc>=A->Iend-A->Istart) {
310: iloc = -1;
311: break;
312: }
313: rcol = A->rcol[iloc];
314: nz = A->nzr[iloc];
315: }
316: if (iloc>=0) {
317: k2 += k1 + 1;
318: rval = A->rval[iloc];
319: *i = A->Istart+iloc;
320: if (!found1 || !found2) {
321: /* Copy row to a new one */
322: rcol2 = rcol;
323: rval2 = rval;
324: PetscCall(PetscMalloc2(nz+1,&rcol,nz+1,&rval));
325: if (found1) {
326: kx = k2;
327: jx = iter->j2;
328: } else {
329: kx = k1;
330: jx = iter->j1;
331: }
332: for (k=0; k<kx; k++) {
333: rcol[k] = rcol2[k];
334: rval[k] = rval2[k];
335: }
336: rcol[kx] = jx;
337: rval[kx] = 0.0;
338: for (k=kx+1; k<nz+1; k++) {
339: rcol[k] = rcol2[k-1];
340: rval[k] = rval2[k-1];
341: }
342: PetscCall(XMatChangeRow(A,*i,nz+1,rcol,rval));
343: }
344: *paij1 = &rval[k1];
345: *paij2 = &rval[k2];
346: iloc++;
347: if (iloc>=A->Iend-A->Istart) iloc = -1;
348: }
349: iter->iloc = iloc;
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode SendrecvRow(XMat A,PetscInt nc1,const PetscInt vj1[],const PetscScalar va1[],PetscInt i2,PetscInt *nc2,PetscInt vj2[],PetscScalar va2[])
355: {
356: PetscInt *ranges,N=A->N;
357: PetscMPIInt rank,naux,len1,len2;
358: MPI_Status st;
360: PetscFunctionBeginUser;
361: PetscCall(XMatGetOwnershipRanges(A,&ranges));
362: PetscAssertAbort(i2>=0 && i2<A->M,A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Incorrect index of row");
363: /* Find owner of row i2 */
364: rank=0;
365: while (ranges[rank+1]<=i2) rank++;
366: /* Send row i1, receive row i2 */
367: PetscCall(PetscMPIIntCast(nc1,&len1));
368: PetscCall(PetscMPIIntCast(N,&len2));
369: PetscCallMPI(MPI_Sendrecv(vj1,len1,MPIU_INT,rank,0,vj2,len2,MPIU_INT,rank,0,A->comm,&st));
370: PetscCallMPI(MPI_Sendrecv(va1,len1,MPIU_SCALAR,rank,0,va2,len2,MPIU_SCALAR,rank,0,A->comm,MPI_STATUS_IGNORE));
371: PetscCallMPI(MPI_Get_count(&st,MPIU_INT,&naux));
372: *nc2 = (PetscInt)naux;
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: static PetscErrorCode PadRows(PetscInt nc1,const PetscInt vj1[],const PetscScalar va1[],PetscInt nc2,const PetscInt vj2[],const PetscScalar va2[],PetscInt *nc,PetscInt **vjj1,PetscScalar **vaa1,PetscInt **vjj2,PetscScalar **vaa2,PetscInt *iwork,PetscScalar *swork)
377: {
378: PetscInt k1,k2,k,n1=nc1+nc2,*vjjaux=iwork;
379: PetscScalar *vaa1aux=swork,*vaa2aux=swork+n1;
381: PetscFunctionBeginUser;
382: *nc=0;
383: for (k1=k2=k=0; k1<nc1 && k2<nc2; ) {
384: if (vj1[k1]<vj2[k2]) {
385: /* Take next element from first row */
386: vaa1aux[k] = va1[k1];
387: vaa2aux[k] = 0.0;
388: vjjaux[k++] = vj1[k1++];
389: } else if (vj1[k1]>vj2[k2]) {
390: /* Take next element from second row */
391: vaa1aux[k] = 0.0;
392: vaa2aux[k] = va2[k2];
393: vjjaux[k++] = vj2[k2++];
394: } else {
395: /* Take next element from both rows */
396: vaa1aux[k] = va1[k1];
397: vaa2aux[k] = va2[k2++];
398: vjjaux[k++] = vj1[k1++];
399: }
400: }
402: /* We reached the end of one of the rows. Continue with the other one */
403: while (k1<nc1) {
404: /* Take next element from first row */
405: vaa1aux[k] = va1[k1];
406: vaa2aux[k] = 0.0;
407: vjjaux[k++] = vj1[k1++];
408: }
409: while (k2<nc2) {
410: /* Take next element from second row */
411: vaa1aux[k] = 0.0;
412: vaa2aux[k] = va2[k2];
413: vjjaux[k++] = vj2[k2++];
414: }
415: *nc=k;
417: /* Copy rows (each row must be allocated separately)*/
418: PetscCall(PetscMalloc2(k,vjj1,k,vaa1));
419: PetscCall(PetscMalloc2(k,vjj2,k,vaa2));
420: PetscCall(PetscArraycpy(*vjj1,vjjaux,k));
421: PetscCall(PetscArraycpy(*vjj2,vjjaux,k));
422: PetscCall(PetscArraycpy(*vaa1,vaa1aux,k));
423: PetscCall(PetscArraycpy(*vaa2,vaa2aux,k));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: static inline void MyAxpby(PetscInt n,PetscScalar a,PetscScalar x[],PetscScalar b,const PetscScalar y[])
428: {
429: PetscInt i;
431: for (i=0; i<n; i++) x[i] = a*x[i]+b*y[i];
432: }
434: /* Apply plane rotation on rows i1, i2 of A */
435: static PetscErrorCode RotateRows(XMat A,PetscInt i1,PetscInt i2,PetscReal c,PetscReal s,PetscInt *nzloc,PetscInt *iwork,PetscScalar *swork)
436: {
437: PetscInt M,N,Istart,Iend,nc,nc1,nc2,*vj1=NULL,*vj2=NULL,*vjj1=NULL,*vjj2=NULL,*iworkx=iwork;
438: PetscBLASInt nc_,one=1;
439: PetscScalar *va1=NULL,*va2=NULL,*vaa1=NULL,*vaa2=NULL,*sworkx=swork;
440: PetscBool i1mine, i2mine;
441: #ifdef TIMING
442: PetscReal t,t0;
443: #endif
445: PetscFunctionBeginUser;
446: PetscCall(XMatGetOwnershipRange(A,&Istart,&Iend));
447: i1mine = (Istart<=i1 && i1<Iend)? PETSC_TRUE: PETSC_FALSE;
448: i2mine = (Istart<=i2 && i2<Iend)? PETSC_TRUE: PETSC_FALSE;
450: if (i1mine || i2mine) {
451: PetscCall(XMatGetSize(A,&M,&N));
453: /* Get local row(s) (at least 1, possibly 2) */
454: time_now(t0);
455: if (i1mine) PetscCall(XMatGetRow(A,i1,&nc1,&vj1,&va1));
456: if (i2mine) PetscCall(XMatGetRow(A,i2,&nc2,&vj2,&va2));
457: time_diff(tt_getr,t0,t,t0);
458: /* Get remote row (at most 1, possibly none) */
459: if (!i2mine) {
460: vj2 = iworkx;
461: iworkx += N;
462: va2 = sworkx;
463: sworkx += N;
464: PetscCall(SendrecvRow(A,nc1,vj1,va1,i2,&nc2,vj2,va2));
465: } else if (!i1mine) {
466: vj1 = iworkx;
467: iworkx += N;
468: va1 = sworkx;
469: sworkx += N;
470: PetscCall(SendrecvRow(A,nc2,vj2,va2,i1,&nc1,vj1,va1));
471: }
473: PetscCall(PadRows(nc1,vj1,va1,nc2,vj2,va2,&nc,&vjj1,&vaa1,&vjj2,&vaa2,iworkx,sworkx));
474: if (i1mine) {
475: *nzloc += nc-nc1;
476: PetscCall(XMatRestoreRow(A,i1,&nc1,&vj1,&va1));
477: }
478: if (i2mine) {
479: *nzloc += nc-nc2;
480: PetscCall(XMatRestoreRow(A,i2,&nc2,&vj2,&va2));
481: }
483: /* Compute rotation and update matrix */
484: if (nc) {
485: PetscCall(PetscBLASIntCast(nc,&nc_));
486: if (i1mine && i2mine) PetscCallBLAS("BLASrot",BLASMIXEDrot_(&nc_,vaa1,&one,vaa2,&one,&c,&s));
487: else if (i1mine) {
488: MyAxpby(nc,c,vaa1,s,vaa2);
489: PetscCall(PetscFree2(vjj2,vaa2));
490: } else {
491: MyAxpby(nc,c,vaa2,-s,vaa1);
492: PetscCall(PetscFree2(vjj1,vaa1));
493: }
494: time_now(t0);
495: if (i1mine) PetscCall(XMatChangeRow(A,i1,nc,vjj1,vaa1));
496: if (i2mine) PetscCall(XMatChangeRow(A,i2,nc,vjj2,vaa2));
497: time_diff(tt_setvalr,t0,t,t0);
498: }
499: }
500: PetscFunctionReturn(PETSC_SUCCESS);
501: }
503: /* Apply plane rotation on columns j1, j2 of A */
504: static PetscErrorCode RotateCols(XMat A,PetscInt j1,PetscInt j2,PetscReal c,PetscReal s,PetscInt *nzloc)
505: {
506: PetscInt i;
507: PetscScalar aux1,aux2,*paij1,*paij2;
508: ColsNzIter iter;
510: PetscFunctionBeginUser;
511: PetscCall(XMatCreateColsNzIter(A,j1,j2,&iter));
512: while (PETSC_TRUE) {
513: PetscCall(ColitNextRow(iter,&i,&paij1,&paij2));
514: if (i<0) break;
515: if (*paij1**paij2==0.0) (*nzloc)++;
516: aux1 = *paij1*c+*paij2*s;
517: aux2 = -*paij1*s+*paij2*c;
518: *paij1 = aux1;
519: *paij2 = aux2;
520: }
521: PetscCall(ColitDestroy(&iter));
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /* Generate a pair of indices i1, i2, where 0 <= i1 < i2 < n */
526: static PetscErrorCode GetIndexPair(PetscRandom rand,PetscInt n,PetscInt *i1,PetscInt *i2)
527: {
528: PetscInt aux;
529: PetscReal x;
531: PetscFunctionBeginUser;
532: PetscCall(PetscRandomGetValueReal(rand,&x));
533: *i1 = (PetscInt)(x*n);
534: PetscCall(PetscRandomGetValueReal(rand,&x));
535: *i2 = (PetscInt)(x*(n-1));
536: if (*i2>=*i1) (*i2)++;
537: if (*i1>*i2) {
538: aux = *i1;
539: *i1 = *i2;
540: *i2 = aux;
541: }
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: int main(int argc,char **argv)
546: {
547: Mat A,Omega; /* operator matrix, signature matrix */
548: XMat XA;
549: SVD svd; /* singular value problem solver context */
550: Vec u,v,vomega,*U,*V;
551: MatType Atype;
552: PetscReal tol,lev1=0.0,lev2=0.0,k=1e3,dens=0.1,log0,loginc,sq,x,angle,c,s;
553: PetscScalar *swork;
554: PetscInt P=110,Q=80,N=100,i,j,Istart,Iend,nconv,nsv,i1,i2,nroth,nrotv,*iwork,nnzwanted,nz,nzloc;
555: PetscBool flg,flgp,flgq,flgn,terse,skiporth=PETSC_FALSE,progress=PETSC_FALSE;
556: PetscRandom rand;
557: MatInfo Ainfo;
558: #ifdef TIMING
559: PetscReal t,t0;
560: #endif
562: PetscFunctionBeginUser;
563: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
565: PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&P,&flgp));
566: PetscCall(PetscOptionsGetInt(NULL,NULL,"-q",&Q,&flgq));
567: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,&flgn));
568: PetscCall(PetscOptionsGetReal(NULL,NULL,"-k",&k,&flg));
569: PetscCall(PetscOptionsGetReal(NULL,NULL,"-d",&dens,&flg));
570: PetscCall(PetscOptionsGetBool(NULL,NULL,"-prog",&progress,&flg));
571: if (!flgp && flgn) P=N;
572: if (!flgq && flgn) Q=N;
573: PetscCheck(P>=N,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p must be >= n");
574: PetscCheck(k>=1.0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter k must be >= 1.0");
575: PetscCheck(dens<=1.0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter d must be <= 1.0");
577: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
578: Build matrix that defines the hyperbolic singular value problem
579: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
581: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHyperbolic singular value problem.\n\n"));
582: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Matrix dimensions: (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT,P,Q,N));
583: PetscCall(PetscPrintf(PETSC_COMM_WORLD,", using signature Omega=blkdiag(I_%" PetscInt_FMT ",-I_%" PetscInt_FMT ")\n\n",P,Q));
585: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
586: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,P+Q,N));
587: PetscCall(MatSetFromOptions(A));
588: PetscCall(MatSetOption(A,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE));
590: /* Set diagonals */
591: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
593: PetscCall(XMatCreate(MPI_COMM_WORLD,&XA,Istart,Iend,P+Q,N));
594: log0 = PetscLogReal(1/k);
595: loginc = -log0/(N-1);
596: sq = PetscSqrtReal(5.0/4.0);
597: for (i=Istart;i<Iend;i++) {
598: if (i<P && i<N) {
599: if (i<Q) x = sq*PetscExpReal(log0+i*loginc);
600: else x = PetscExpReal(log0+i*loginc);
601: PetscCall(XMatSetValue(XA,i,i,x,INSERT_VALUES));
602: } else if (i>=P && i-P<N) {
603: j = i-P;
604: PetscCall(XMatSetValue(XA,i,j,PetscExpReal(log0+j*loginc),INSERT_VALUES));
605: }
606: }
608: /* Apply plane rotations */
609: nnzwanted = dens*(P+Q)*N;
610: PetscCall(XMatGetNumberNonzeros(XA,&nzloc));
611: PetscCallMPI(MPI_Allreduce(&nzloc,&nz,1,MPIU_INT,MPI_SUM,MPI_COMM_WORLD));
612: PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand));
613: PetscCall(PetscRandomSetFromOptions(rand));
614: nroth = nrotv = 0;
615: PetscCall(PetscMalloc2(2*N,&iwork,3*N,&swork));
616: time_now(t0);
617: while (nz<0.95*nnzwanted) {
618: /* Plane rotation on 2 of the top P rows*/
619: PetscCall(PetscRandomGetValueReal(rand,&angle));
620: c=PetscCosReal(angle);
621: s=PetscSinReal(angle);
622: PetscCall(GetIndexPair(rand,P,&i1,&i2));
623: PetscCall(RotateRows(XA,i1,i2,c,s,&nzloc,iwork,swork));
624: time_diff(tt_rotr,t0,t,t0);
625: nroth++;
626: /* Plane rotation on 2 of the bottom Q rows*/
627: if (Q>1) {
628: PetscCall(PetscRandomGetValueReal(rand,&angle));
629: c=PetscCosReal(angle);
630: s=PetscSinReal(angle);
631: PetscCall(GetIndexPair(rand,Q,&i1,&i2));
632: PetscCall(RotateRows(XA,P+i1,P+i2,c,s,&nzloc,iwork,swork));
633: time_diff(tt_rotr,t0,t,t0);
634: nroth++;
635: }
636: /* Plane rotation on 2 columns */
637: PetscCall(PetscRandomGetValueReal(rand,&angle));
638: c=PetscCosReal(angle);
639: s=PetscSinReal(angle);
640: PetscCall(GetIndexPair(rand,N,&i1,&i2));
641: PetscCall(RotateCols(XA,i1,i2,c,s,&nzloc));
642: time_diff(tt_rotc,t0,t,t0);
643: nrotv++;
644: /* Update total number of non-zeros */
645: PetscCallMPI(MPI_Allreduce(&nzloc,&nz,1,MPIU_INT,MPI_SUM,MPI_COMM_WORLD));
646: if (progress) PetscCall(PetscFPrintf(MPI_COMM_WORLD,stderr,"\rProgress: %.2f%%",(double)nz/nnzwanted*100));
647: }
648: if (progress) PetscCall(PetscFPrintf(MPI_COMM_WORLD,stderr,"\r"));
650: PetscCall(PetscFree2(iwork,swork));
651: PetscCall(PetscRandomDestroy(&rand));
653: PetscCall(XMatConvert(XA,A));
654: PetscCall(XMatDestroy(&XA));
655: PetscCall(MatGetInfo(A,MAT_GLOBAL_SUM,&Ainfo));
656: PetscCall(PetscPrintf(MPI_COMM_WORLD," Matrix created. %" PetscInt_FMT " rotations applied. nnz=%.0f. Density: %.4g\n\n",nroth+nrotv,Ainfo.nz_used,(double)Ainfo.nz_used/(1.0*(P+Q)*N)));
658: #ifdef TIMING
659: PetscCall(PetscPrintf(MPI_COMM_WORLD,"#T rot-rows: %6.3f get-rows: %6.3f setval-rows: %6.3f assm-rows: %.3f\n",tt_rotr,tt_getr,tt_setvalr,tt_assmr));
660: PetscCall(PetscPrintf(MPI_COMM_WORLD,"#T rot cols: %6.3f get-cols: %6.3f setval-cols: %6.3f assm-cols: %.3f\n",tt_rotc,tt_getc,tt_setvalc,tt_assmc));
661: #endif
663: /* scale by 0.5 the lower Q rows of A */
664: PetscCall(MatCreateVecs(A,NULL,&vomega));
665: PetscCall(VecSet(vomega,1.0));
666: for (i=PetscMax(P,Istart);i<Iend;i++) {
667: PetscCall(VecSetValue(vomega,i,0.5,INSERT_VALUES));
668: }
669: PetscCall(VecAssemblyBegin(vomega));
670: PetscCall(VecAssemblyEnd(vomega));
671: PetscCall(MatDiagonalScale(A,vomega,NULL));
673: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
674: Create the signature
675: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
677: PetscCall(VecSet(vomega,1.0));
678: PetscCall(VecGetOwnershipRange(vomega,&Istart,&Iend));
679: for (i=Istart;i<Iend;i++) {
680: if (i>=P) PetscCall(VecSetValue(vomega,i,-1.0,INSERT_VALUES));
681: }
682: PetscCall(VecAssemblyBegin(vomega));
683: PetscCall(VecAssemblyEnd(vomega));
685: PetscCall(MatGetType(A,&Atype));
686: PetscCall(MatCreate(PETSC_COMM_WORLD,&Omega));
687: PetscCall(MatSetSizes(Omega,PETSC_DECIDE,PETSC_DECIDE,P+Q,P+Q));
688: PetscCall(MatSetType(Omega,Atype));
689: PetscCall(MatDiagonalSet(Omega,vomega,INSERT_VALUES));
691: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
692: Create the singular value solver and set various options
693: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
695: PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
697: PetscCall(SVDSetOperators(svd,A,NULL));
698: PetscCall(SVDSetSignature(svd,vomega));
699: PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
701: PetscCall(SVDSetFromOptions(svd));
703: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
704: Solve the problem, display solution
705: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
707: PetscCall(MatCreateVecs(A,&v,&u));
708: PetscCall(SVDSolve(svd));
710: /* show detailed info unless -terse option is given by user */
711: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
712: if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
713: else {
714: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
715: PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
716: PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
717: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
718: }
720: /* check orthogonality */
721: PetscCall(PetscOptionsGetBool(NULL,NULL,"-skiporth",&skiporth,NULL));
722: PetscCall(SVDGetConverged(svd,&nconv));
723: PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
724: nconv = PetscMin(nconv,nsv);
725: if (nconv>0 && !skiporth) {
726: PetscCall(SVDGetTolerances(svd,&tol,NULL));
727: PetscCall(VecDuplicateVecs(u,nconv,&U));
728: PetscCall(VecDuplicateVecs(v,nconv,&V));
729: for (i=0;i<nconv;i++) PetscCall(SVDGetSingularTriplet(svd,i,NULL,U[i],V[i]));
730: PetscCall(VecCheckOrthonormality(U,nconv,NULL,nconv,Omega,NULL,&lev1));
731: PetscCall(VecCheckOrthonormality(V,nconv,NULL,nconv,NULL,NULL,&lev2));
732: if (lev1+lev2<20*tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n"));
733: else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g (U) %g (V)\n",(double)lev1,(double)lev2));
734: PetscCall(VecDestroyVecs(nconv,&U));
735: PetscCall(VecDestroyVecs(nconv,&V));
736: }
737: PetscCall(VecDestroy(&u));
738: PetscCall(VecDestroy(&v));
740: /* free work space */
741: PetscCall(SVDDestroy(&svd));
742: PetscCall(MatDestroy(&Omega));
743: PetscCall(MatDestroy(&A));
744: PetscCall(VecDestroy(&vomega));
745: PetscCall(SlepcFinalize());
746: return 0;
747: }
749: /*TEST
751: testset:
752: args: -svd_nsv 5 -d .15 -terse
753: output_file: output/ex53_1.out
754: nsize: {{1 2}}
755: test:
756: args: -svd_type trlanczos -ds_parallel {{redundant synchronized}}
757: suffix: 1_trlanczos
758: test:
759: args: -svd_type cross
760: suffix: 1_cross
761: test:
762: args: -svd_type cyclic
763: requires: !single
764: suffix: 1_cyclic
766: TEST*/