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[] = "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";
19 :
20 : #include <slepcsvd.h>
21 : #include <slepcblaslapack.h>
22 :
23 : /* The problem matrix is the (p+q)-by-n matrix
24 :
25 : A = | Q1 * D1 * U |
26 : | 1/2*Q2 * D2 * U |
27 :
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.
35 :
36 : The signature matrix is Omega = diag(I_p,-I_q)
37 :
38 : The hyperbolic singular values of A are equal to the elements of vector d.
39 :
40 : A is generated by random plane rotations applied to the diagonal matrices
41 : D1 and D2
42 :
43 : This test case is based on the paper:
44 :
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 : */
49 :
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
59 :
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 : };
69 :
70 : typedef struct _p_XMat *XMat;
71 :
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 : };
78 :
79 : typedef struct _p_ColsNzIter *ColsNzIter;
80 :
81 12 : static PetscErrorCode XMatCreate(MPI_Comm comm,XMat *A,PetscInt Istart,PetscInt Iend,PetscInt M,PetscInt N)
82 : {
83 12 : PetscInt i;
84 12 : PetscMPIInt np;
85 :
86 12 : PetscFunctionBeginUser;
87 12 : PetscCall(PetscNew(A));
88 12 : (*A)->M = M;
89 12 : (*A)->N = N;
90 12 : (*A)->Istart = Istart;
91 12 : (*A)->Iend = Iend;
92 12 : (*A)->comm = comm;
93 12 : PetscCallMPI(MPI_Comm_size(comm,&np));
94 12 : PetscCall(PetscMalloc4(Iend-Istart,&(*A)->nzr,Iend-Istart,&(*A)->rcol,Iend-Istart,&(*A)->rval,np+1,&(*A)->ranges));
95 1532 : for (i=0; i<Iend-Istart; i++) {
96 1520 : (*A)->nzr[i] = 0;
97 1520 : (*A)->rcol[i] = NULL;
98 1520 : (*A)->rval[i] = NULL;
99 : }
100 24 : PetscCallMPI(MPI_Allgather(&Istart,1,MPIU_INT,(*A)->ranges,1,MPIU_INT,comm));
101 12 : (*A)->ranges[np]=M;
102 12 : PetscFunctionReturn(PETSC_SUCCESS);
103 : }
104 :
105 12 : static PetscErrorCode XMatDestroy(XMat *A)
106 : {
107 12 : PetscInt m,i;
108 :
109 12 : PetscFunctionBeginUser;
110 12 : m = (*A)->Iend-(*A)->Istart;
111 1532 : for (i=0; i<m; i++) {
112 1520 : if ((*A)->nzr[i]>0) {
113 1520 : PetscCall(PetscFree2((*A)->rcol[i],(*A)->rval[i]));
114 : }
115 : }
116 12 : PetscCall(PetscFree4((*A)->nzr,(*A)->rcol,(*A)->rval,(*A)->ranges));
117 12 : PetscCall(PetscFree(*A));
118 12 : *A = NULL;
119 12 : PetscFunctionReturn(PETSC_SUCCESS);
120 : }
121 :
122 1264 : static PetscErrorCode XMatGetSize(XMat A,PetscInt *M,PetscInt *N)
123 : {
124 1264 : PetscFunctionBeginUser;
125 1264 : *M = A->M;
126 1264 : *N = A->N;
127 1264 : PetscFunctionReturn(PETSC_SUCCESS);
128 : }
129 :
130 1776 : static PetscErrorCode XMatGetOwnershipRange(XMat A,PetscInt *Istart,PetscInt *Iend)
131 : {
132 1776 : PetscFunctionBeginUser;
133 1776 : *Istart = A->Istart;
134 1776 : *Iend = A->Iend;
135 1776 : PetscFunctionReturn(PETSC_SUCCESS);
136 : }
137 :
138 160 : static PetscErrorCode XMatGetOwnershipRanges(XMat A,PetscInt **ranges)
139 : {
140 160 : PetscFunctionBeginUser;
141 160 : *ranges = A->ranges;
142 160 : PetscFunctionReturn(PETSC_SUCCESS);
143 : }
144 :
145 2368 : static PetscErrorCode XMatGetRow(XMat A,PetscInt i,PetscInt *nc,PetscInt **vj,PetscScalar **va)
146 : {
147 2368 : PetscInt iloc=i-A->Istart;
148 :
149 2368 : PetscFunctionBeginUser;
150 2368 : *nc = A->nzr[iloc];
151 2368 : *vj = A->rcol[iloc];
152 2368 : *va = A->rval[iloc];
153 2368 : PetscFunctionReturn(PETSC_SUCCESS);
154 : }
155 :
156 2368 : static PetscErrorCode XMatRestoreRow(XMat A,PetscInt i,PetscInt *nc,PetscInt **vj,PetscScalar **va)
157 : {
158 2368 : PetscFunctionBeginUser;
159 2368 : *nc = 0;
160 2368 : *vj = NULL;
161 2368 : *va = NULL;
162 2368 : PetscFunctionReturn(PETSC_SUCCESS);
163 : }
164 :
165 12160 : static PetscErrorCode XMatChangeRow(XMat A,PetscInt i,PetscInt nc,PetscInt vj[],PetscScalar va[])
166 : {
167 12160 : PetscInt iloc = i-A->Istart;
168 :
169 12160 : PetscFunctionBeginUser;
170 12160 : if (A->nzr[iloc]>0) {
171 12088 : PetscCall(PetscFree2(A->rcol[iloc],A->rval[iloc]));
172 : }
173 12160 : A->nzr[iloc] = nc;
174 12160 : A->rcol[iloc] = vj;
175 12160 : A->rval[iloc] = va;
176 12160 : PetscFunctionReturn(PETSC_SUCCESS);
177 : }
178 :
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 226400 : static PetscBool getcolpos(const PetscInt rcol[],PetscInt nz,PetscInt j,PetscInt *k)
183 : {
184 226400 : PetscInt ka,kb;
185 226400 : if (nz==0) {
186 23952 : *k = 0;
187 23952 : return PETSC_FALSE;
188 : }
189 202448 : if (j<=rcol[0]) {
190 67504 : *k = 0;
191 67504 : return (rcol[0]==j)? PETSC_TRUE: PETSC_FALSE;
192 : }
193 134944 : if (rcol[nz-1]<=j) {
194 45680 : *k = nz-(rcol[nz-1]==j);
195 45680 : return (rcol[nz-1]==j)? PETSC_TRUE: PETSC_FALSE;
196 : }
197 : /* Bisection: rcol[ka] < j < rcol[kb] */
198 89264 : ka = 0; kb = nz-1;
199 277048 : while (ka+1<kb) {
200 195048 : *k = (ka+kb)/2;
201 195048 : if (rcol[*k]<=j) {
202 98032 : if (rcol[*k]==j) return PETSC_TRUE;
203 : else ka = *k;
204 : } else kb = *k;
205 : }
206 82000 : *k = kb;
207 82000 : return PETSC_FALSE;
208 : }
209 :
210 : /* Only local elements can be set */
211 1440 : static PetscErrorCode XMatSetValue(XMat A,PetscInt i,PetscInt j,PetscScalar x,InsertMode dummy)
212 : {
213 1440 : PetscInt nz,iloc,k,kx,*col1,*col2;
214 1440 : PetscScalar *val1,*val2;
215 :
216 1440 : PetscFunctionBeginUser;
217 1440 : PetscCheck(i>=A->Istart && i<A->Iend && j>=0 && j<A->N,A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Incorrect index of row/column");
218 1440 : iloc = i-A->Istart;
219 1440 : nz = A->nzr[iloc];
220 1440 : col1 = A->rcol[iloc];
221 1440 : val1 = A->rval[iloc];
222 1440 : if (!getcolpos(col1,nz,j,&k)) {
223 1440 : A->nzr[iloc]++;
224 : /* Replace row */
225 1440 : PetscCall(PetscMalloc2(nz+1,&col2,nz+1,&val2));
226 1440 : A->rcol[iloc] = col2;
227 1440 : A->rval[iloc] = val2;
228 1440 : for (kx=0; kx<k; kx++) {
229 0 : col2[kx] = col1[kx];
230 0 : val2[kx] = val1[kx];
231 : }
232 1440 : for (; kx<nz; kx++) {
233 0 : col2[kx+1] = col1[kx];
234 0 : val2[kx+1] = val1[kx];
235 : }
236 1440 : if (nz>0) {
237 0 : PetscCall(PetscFree2(col1,val1));
238 : }
239 1440 : col1 = col2;
240 1440 : val1 = val2;
241 : }
242 1440 : col1[k] = j;
243 1440 : val1[k] = x;
244 1440 : PetscFunctionReturn(PETSC_SUCCESS);
245 : }
246 :
247 : /* Convert to PETSc Mat */
248 12 : static PetscErrorCode XMatConvert(XMat XA,Mat A)
249 : {
250 12 : PetscInt i,iloc;
251 :
252 12 : PetscFunctionBeginUser;
253 1532 : for (i=XA->Istart; i<XA->Iend; i++) {
254 1520 : iloc = i-XA->Istart;
255 1520 : PetscCall(MatSetValues(A,1,&i,XA->nzr[iloc],XA->rcol[iloc],XA->rval[iloc],INSERT_VALUES));
256 : }
257 12 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
258 12 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
259 12 : PetscFunctionReturn(PETSC_SUCCESS);
260 : }
261 :
262 : /* gets the number of nonzeros in the matrix on this MPI rank */
263 12 : static PetscErrorCode XMatGetNumberNonzeros(XMat A,PetscInt *nz)
264 : {
265 12 : PetscInt i;
266 :
267 12 : PetscFunctionBeginUser;
268 12 : *nz = 0;
269 1532 : for (i=0; i<A->Iend-A->Istart; i++) *nz += A->nzr[i];
270 12 : PetscFunctionReturn(PETSC_SUCCESS);
271 : }
272 :
273 888 : static PetscErrorCode XMatCreateColsNzIter(XMat A,PetscInt j1,PetscInt j2,ColsNzIter *iter)
274 : {
275 888 : PetscFunctionBeginUser;
276 888 : PetscCall(PetscNew(iter));
277 888 : (*iter)->A = A;
278 888 : (*iter)->j1 = j1;
279 888 : (*iter)->j2 = j2;
280 888 : (*iter)->iloc = 0;
281 888 : PetscFunctionReturn(PETSC_SUCCESS);
282 : }
283 :
284 888 : static PetscErrorCode ColitDestroy(ColsNzIter *iter)
285 : {
286 888 : PetscFunctionBeginUser;
287 888 : PetscCall(PetscFree(*iter));
288 888 : PetscFunctionReturn(PETSC_SUCCESS);
289 : }
290 :
291 11560 : static PetscErrorCode ColitNextRow(ColsNzIter iter,PetscInt *i,PetscScalar **paij1,PetscScalar **paij2)
292 : {
293 11560 : PetscInt iloc = iter->iloc,*rcol,nz,k,k1,k2,kx,jx,*rcol2;
294 11560 : PetscScalar *rval,*rval2;
295 11560 : PetscBool found1, found2;
296 11560 : XMat A=iter->A;
297 :
298 11560 : PetscFunctionBeginUser;
299 11560 : *i = -1;
300 11560 : *paij1 = *paij2 = NULL;
301 11560 : if (iloc>=0) {
302 11488 : rcol = A->rcol[iloc];
303 11488 : nz = A->nzr[iloc];
304 213472 : while (PETSC_TRUE) {
305 112480 : found1 = getcolpos(rcol,nz,iter->j1,&k1);
306 112480 : found2 = getcolpos(PetscSafePointerPlusOffset(rcol,k1+found1),nz-k1-found1,iter->j2,&k2);
307 112480 : if (found1 || found2) break;
308 101808 : iloc++;
309 101808 : if (iloc>=A->Iend-A->Istart) {
310 : iloc = -1;
311 : break;
312 : }
313 100992 : rcol = A->rcol[iloc];
314 100992 : nz = A->nzr[iloc];
315 : }
316 11488 : if (iloc>=0) {
317 10672 : k2 += k1 + 1;
318 10672 : rval = A->rval[iloc];
319 10672 : *i = A->Istart+iloc;
320 10672 : if (!found1 || !found2) {
321 : /* Copy row to a new one */
322 9808 : rcol2 = rcol;
323 9808 : rval2 = rval;
324 9808 : PetscCall(PetscMalloc2(nz+1,&rcol,nz+1,&rval));
325 9808 : if (found1) {
326 5472 : kx = k2;
327 5472 : jx = iter->j2;
328 : } else {
329 4336 : kx = k1;
330 4336 : jx = iter->j1;
331 : }
332 63976 : for (k=0; k<kx; k++) {
333 54168 : rcol[k] = rcol2[k];
334 54168 : rval[k] = rval2[k];
335 : }
336 9808 : rcol[kx] = jx;
337 9808 : rval[kx] = 0.0;
338 72784 : for (k=kx+1; k<nz+1; k++) {
339 62976 : rcol[k] = rcol2[k-1];
340 62976 : rval[k] = rval2[k-1];
341 : }
342 9808 : PetscCall(XMatChangeRow(A,*i,nz+1,rcol,rval));
343 : }
344 10672 : *paij1 = &rval[k1];
345 10672 : *paij2 = &rval[k2];
346 10672 : iloc++;
347 10672 : if (iloc>=A->Iend-A->Istart) iloc = -1;
348 : }
349 11488 : iter->iloc = iloc;
350 : }
351 11560 : PetscFunctionReturn(PETSC_SUCCESS);
352 : }
353 :
354 160 : static PetscErrorCode SendrecvRow(XMat A,PetscInt nc1,const PetscInt vj1[],const PetscScalar va1[],PetscInt i2,PetscInt *nc2,PetscInt vj2[],PetscScalar va2[])
355 : {
356 160 : PetscInt *ranges,N=A->N;
357 160 : PetscMPIInt rank,naux,len1,len2;
358 160 : MPI_Status st;
359 :
360 160 : PetscFunctionBeginUser;
361 160 : PetscCall(XMatGetOwnershipRanges(A,&ranges));
362 160 : PetscAssertAbort(i2>=0 && i2<A->M,A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Incorrect index of row");
363 : /* Find owner of row i2 */
364 160 : rank=0;
365 240 : while (ranges[rank+1]<=i2) rank++;
366 : /* Send row i1, receive row i2 */
367 160 : PetscCall(PetscMPIIntCast(nc1,&len1));
368 160 : PetscCall(PetscMPIIntCast(N,&len2));
369 160 : PetscCallMPI(MPI_Sendrecv(vj1,len1,MPIU_INT,rank,0,vj2,len2,MPIU_INT,rank,0,A->comm,&st));
370 160 : PetscCallMPI(MPI_Sendrecv(va1,len1,MPIU_SCALAR,rank,0,va2,len2,MPIU_SCALAR,rank,0,A->comm,MPI_STATUS_IGNORE));
371 160 : PetscCallMPI(MPI_Get_count(&st,MPIU_INT,&naux));
372 160 : *nc2 = (PetscInt)naux;
373 160 : PetscFunctionReturn(PETSC_SUCCESS);
374 : }
375 :
376 1264 : 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 1264 : PetscInt k1,k2,k,n1=nc1+nc2,*vjjaux=iwork;
379 1264 : PetscScalar *vaa1aux=swork,*vaa2aux=swork+n1;
380 :
381 1264 : PetscFunctionBeginUser;
382 1264 : *nc=0;
383 10560 : for (k1=k2=k=0; k1<nc1 && k2<nc2; ) {
384 9296 : if (vj1[k1]<vj2[k2]) {
385 : /* Take next element from first row */
386 4324 : vaa1aux[k] = va1[k1];
387 4324 : vaa2aux[k] = 0.0;
388 4324 : vjjaux[k++] = vj1[k1++];
389 4972 : } else if (vj1[k1]>vj2[k2]) {
390 : /* Take next element from second row */
391 4312 : vaa1aux[k] = 0.0;
392 4312 : vaa2aux[k] = va2[k2];
393 4312 : vjjaux[k++] = vj2[k2++];
394 : } else {
395 : /* Take next element from both rows */
396 660 : vaa1aux[k] = va1[k1];
397 660 : vaa2aux[k] = va2[k2++];
398 660 : vjjaux[k++] = vj1[k1++];
399 : }
400 : }
401 :
402 : /* We reached the end of one of the rows. Continue with the other one */
403 2556 : while (k1<nc1) {
404 : /* Take next element from first row */
405 1292 : vaa1aux[k] = va1[k1];
406 1292 : vaa2aux[k] = 0.0;
407 1292 : vjjaux[k++] = vj1[k1++];
408 : }
409 2580 : while (k2<nc2) {
410 : /* Take next element from second row */
411 1316 : vaa1aux[k] = 0.0;
412 1316 : vaa2aux[k] = va2[k2];
413 1316 : vjjaux[k++] = vj2[k2++];
414 : }
415 1264 : *nc=k;
416 :
417 : /* Copy rows (each row must be allocated separately)*/
418 1264 : PetscCall(PetscMalloc2(k,vjj1,k,vaa1));
419 1264 : PetscCall(PetscMalloc2(k,vjj2,k,vaa2));
420 1264 : PetscCall(PetscArraycpy(*vjj1,vjjaux,k));
421 1264 : PetscCall(PetscArraycpy(*vjj2,vjjaux,k));
422 1264 : PetscCall(PetscArraycpy(*vaa1,vaa1aux,k));
423 1264 : PetscCall(PetscArraycpy(*vaa2,vaa2aux,k));
424 1264 : PetscFunctionReturn(PETSC_SUCCESS);
425 : }
426 :
427 160 : static inline void MyAxpby(PetscInt n,PetscScalar a,PetscScalar x[],PetscScalar b,const PetscScalar y[])
428 : {
429 160 : PetscInt i;
430 :
431 1280 : for (i=0; i<n; i++) x[i] = a*x[i]+b*y[i];
432 : }
433 :
434 : /* Apply plane rotation on rows i1, i2 of A */
435 1776 : static PetscErrorCode RotateRows(XMat A,PetscInt i1,PetscInt i2,PetscReal c,PetscReal s,PetscInt *nzloc,PetscInt *iwork,PetscScalar *swork)
436 : {
437 1776 : PetscInt M,N,Istart,Iend,nc,nc1,nc2,*vj1=NULL,*vj2=NULL,*vjj1=NULL,*vjj2=NULL,*iworkx=iwork;
438 1776 : PetscBLASInt nc_,one=1;
439 1776 : PetscScalar *va1=NULL,*va2=NULL,*vaa1=NULL,*vaa2=NULL,*sworkx=swork;
440 1776 : PetscBool i1mine, i2mine;
441 : #ifdef TIMING
442 : PetscReal t,t0;
443 : #endif
444 :
445 1776 : PetscFunctionBeginUser;
446 1776 : PetscCall(XMatGetOwnershipRange(A,&Istart,&Iend));
447 1776 : i1mine = (Istart<=i1 && i1<Iend)? PETSC_TRUE: PETSC_FALSE;
448 1776 : i2mine = (Istart<=i2 && i2<Iend)? PETSC_TRUE: PETSC_FALSE;
449 :
450 1776 : if (i1mine || i2mine) {
451 1264 : PetscCall(XMatGetSize(A,&M,&N));
452 :
453 : /* Get local row(s) (at least 1, possibly 2) */
454 1264 : time_now(t0);
455 1264 : if (i1mine) PetscCall(XMatGetRow(A,i1,&nc1,&vj1,&va1));
456 1264 : if (i2mine) PetscCall(XMatGetRow(A,i2,&nc2,&vj2,&va2));
457 1264 : time_diff(tt_getr,t0,t,t0);
458 : /* Get remote row (at most 1, possibly none) */
459 1264 : if (!i2mine) {
460 80 : vj2 = iworkx;
461 80 : iworkx += N;
462 80 : va2 = sworkx;
463 80 : sworkx += N;
464 80 : PetscCall(SendrecvRow(A,nc1,vj1,va1,i2,&nc2,vj2,va2));
465 1184 : } else if (!i1mine) {
466 80 : vj1 = iworkx;
467 80 : iworkx += N;
468 80 : va1 = sworkx;
469 80 : sworkx += N;
470 80 : PetscCall(SendrecvRow(A,nc2,vj2,va2,i1,&nc1,vj1,va1));
471 : }
472 :
473 1264 : PetscCall(PadRows(nc1,vj1,va1,nc2,vj2,va2,&nc,&vjj1,&vaa1,&vjj2,&vaa2,iworkx,sworkx));
474 1264 : if (i1mine) {
475 1184 : *nzloc += nc-nc1;
476 1184 : PetscCall(XMatRestoreRow(A,i1,&nc1,&vj1,&va1));
477 : }
478 1264 : if (i2mine) {
479 1184 : *nzloc += nc-nc2;
480 1184 : PetscCall(XMatRestoreRow(A,i2,&nc2,&vj2,&va2));
481 : }
482 :
483 : /* Compute rotation and update matrix */
484 1264 : if (nc) {
485 1256 : PetscCall(PetscBLASIntCast(nc,&nc_));
486 1256 : if (i1mine && i2mine) PetscCallBLAS("BLASrot",BLASMIXEDrot_(&nc_,vaa1,&one,vaa2,&one,&c,&s));
487 160 : else if (i1mine) {
488 80 : MyAxpby(nc,c,vaa1,s,vaa2);
489 80 : PetscCall(PetscFree2(vjj2,vaa2));
490 : } else {
491 80 : MyAxpby(nc,c,vaa2,-s,vaa1);
492 80 : PetscCall(PetscFree2(vjj1,vaa1));
493 : }
494 1256 : time_now(t0);
495 1256 : if (i1mine) PetscCall(XMatChangeRow(A,i1,nc,vjj1,vaa1));
496 1256 : if (i2mine) PetscCall(XMatChangeRow(A,i2,nc,vjj2,vaa2));
497 1776 : time_diff(tt_setvalr,t0,t,t0);
498 : }
499 : }
500 1776 : PetscFunctionReturn(PETSC_SUCCESS);
501 : }
502 :
503 : /* Apply plane rotation on columns j1, j2 of A */
504 888 : static PetscErrorCode RotateCols(XMat A,PetscInt j1,PetscInt j2,PetscReal c,PetscReal s,PetscInt *nzloc)
505 : {
506 888 : PetscInt i;
507 888 : PetscScalar aux1,aux2,*paij1,*paij2;
508 888 : ColsNzIter iter;
509 :
510 888 : PetscFunctionBeginUser;
511 888 : PetscCall(XMatCreateColsNzIter(A,j1,j2,&iter));
512 22232 : while (PETSC_TRUE) {
513 11560 : PetscCall(ColitNextRow(iter,&i,&paij1,&paij2));
514 11560 : if (i<0) break;
515 10672 : if (*paij1**paij2==0.0) (*nzloc)++;
516 10672 : aux1 = *paij1*c+*paij2*s;
517 10672 : aux2 = -*paij1*s+*paij2*c;
518 10672 : *paij1 = aux1;
519 10672 : *paij2 = aux2;
520 : }
521 888 : PetscCall(ColitDestroy(&iter));
522 888 : PetscFunctionReturn(PETSC_SUCCESS);
523 : }
524 :
525 : /* Generate a pair of indices i1, i2, where 0 <= i1 < i2 < n */
526 2664 : static PetscErrorCode GetIndexPair(PetscRandom rand,PetscInt n,PetscInt *i1,PetscInt *i2)
527 : {
528 2664 : PetscInt aux;
529 2664 : PetscReal x;
530 :
531 2664 : PetscFunctionBeginUser;
532 2664 : PetscCall(PetscRandomGetValueReal(rand,&x));
533 2664 : *i1 = (PetscInt)(x*n);
534 2664 : PetscCall(PetscRandomGetValueReal(rand,&x));
535 2664 : *i2 = (PetscInt)(x*(n-1));
536 2664 : if (*i2>=*i1) (*i2)++;
537 2664 : if (*i1>*i2) {
538 1332 : aux = *i1;
539 1332 : *i1 = *i2;
540 1332 : *i2 = aux;
541 : }
542 2664 : PetscFunctionReturn(PETSC_SUCCESS);
543 : }
544 :
545 12 : int main(int argc,char **argv)
546 : {
547 12 : Mat A,Omega; /* operator matrix, signature matrix */
548 12 : XMat XA;
549 12 : SVD svd; /* singular value problem solver context */
550 12 : Vec u,v,vomega,*U,*V;
551 12 : MatType Atype;
552 12 : PetscReal tol,lev1=0.0,lev2=0.0,k=1e3,dens=0.1,log0,loginc,sq,x,angle,c,s;
553 12 : PetscScalar *swork;
554 12 : PetscInt P=110,Q=80,N=100,i,j,Istart,Iend,nconv,nsv,i1,i2,nroth,nrotv,*iwork,nnzwanted,nz,nzloc;
555 12 : PetscBool flg,flgp,flgq,flgn,terse,skiporth=PETSC_FALSE,progress=PETSC_FALSE;
556 12 : PetscRandom rand;
557 12 : MatInfo Ainfo;
558 : #ifdef TIMING
559 : PetscReal t,t0;
560 : #endif
561 :
562 12 : PetscFunctionBeginUser;
563 12 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
564 :
565 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&P,&flgp));
566 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-q",&Q,&flgq));
567 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,&flgn));
568 12 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-k",&k,&flg));
569 12 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-d",&dens,&flg));
570 12 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-prog",&progress,&flg));
571 12 : if (!flgp && flgn) P=N;
572 12 : if (!flgq && flgn) Q=N;
573 12 : PetscCheck(P>=N,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p must be >= n");
574 12 : PetscCheck(k>=1.0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter k must be >= 1.0");
575 12 : PetscCheck(dens<=1.0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter d must be <= 1.0");
576 :
577 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
578 : Build matrix that defines the hyperbolic singular value problem
579 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
580 :
581 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHyperbolic singular value problem.\n\n"));
582 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Matrix dimensions: (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT,P,Q,N));
583 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,", using signature Omega=blkdiag(I_%" PetscInt_FMT ",-I_%" PetscInt_FMT ")\n\n",P,Q));
584 :
585 12 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
586 12 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,P+Q,N));
587 12 : PetscCall(MatSetFromOptions(A));
588 12 : PetscCall(MatSetOption(A,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE));
589 :
590 : /* Set diagonals */
591 12 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
592 :
593 12 : PetscCall(XMatCreate(MPI_COMM_WORLD,&XA,Istart,Iend,P+Q,N));
594 12 : log0 = PetscLogReal(1/k);
595 12 : loginc = -log0/(N-1);
596 12 : sq = PetscSqrtReal(5.0/4.0);
597 1532 : for (i=Istart;i<Iend;i++) {
598 1520 : if (i<P && i<N) {
599 800 : if (i<Q) x = sq*PetscExpReal(log0+i*loginc);
600 160 : else x = PetscExpReal(log0+i*loginc);
601 800 : PetscCall(XMatSetValue(XA,i,i,x,INSERT_VALUES));
602 720 : } else if (i>=P && i-P<N) {
603 640 : j = i-P;
604 1520 : PetscCall(XMatSetValue(XA,i,j,PetscExpReal(log0+j*loginc),INSERT_VALUES));
605 : }
606 : }
607 :
608 : /* Apply plane rotations */
609 12 : nnzwanted = dens*(P+Q)*N;
610 12 : PetscCall(XMatGetNumberNonzeros(XA,&nzloc));
611 24 : PetscCallMPI(MPI_Allreduce(&nzloc,&nz,1,MPIU_INT,MPI_SUM,MPI_COMM_WORLD));
612 12 : PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand));
613 12 : PetscCall(PetscRandomSetFromOptions(rand));
614 12 : nroth = nrotv = 0;
615 12 : PetscCall(PetscMalloc2(2*N,&iwork,3*N,&swork));
616 : time_now(t0);
617 900 : while (nz<0.95*nnzwanted) {
618 : /* Plane rotation on 2 of the top P rows*/
619 888 : PetscCall(PetscRandomGetValueReal(rand,&angle));
620 888 : c=PetscCosReal(angle);
621 888 : s=PetscSinReal(angle);
622 888 : PetscCall(GetIndexPair(rand,P,&i1,&i2));
623 888 : PetscCall(RotateRows(XA,i1,i2,c,s,&nzloc,iwork,swork));
624 888 : time_diff(tt_rotr,t0,t,t0);
625 888 : nroth++;
626 : /* Plane rotation on 2 of the bottom Q rows*/
627 888 : if (Q>1) {
628 888 : PetscCall(PetscRandomGetValueReal(rand,&angle));
629 888 : c=PetscCosReal(angle);
630 888 : s=PetscSinReal(angle);
631 888 : PetscCall(GetIndexPair(rand,Q,&i1,&i2));
632 888 : PetscCall(RotateRows(XA,P+i1,P+i2,c,s,&nzloc,iwork,swork));
633 888 : time_diff(tt_rotr,t0,t,t0);
634 888 : nroth++;
635 : }
636 : /* Plane rotation on 2 columns */
637 888 : PetscCall(PetscRandomGetValueReal(rand,&angle));
638 888 : c=PetscCosReal(angle);
639 888 : s=PetscSinReal(angle);
640 888 : PetscCall(GetIndexPair(rand,N,&i1,&i2));
641 888 : PetscCall(RotateCols(XA,i1,i2,c,s,&nzloc));
642 888 : time_diff(tt_rotc,t0,t,t0);
643 888 : nrotv++;
644 : /* Update total number of non-zeros */
645 1776 : PetscCallMPI(MPI_Allreduce(&nzloc,&nz,1,MPIU_INT,MPI_SUM,MPI_COMM_WORLD));
646 1788 : if (progress) PetscCall(PetscFPrintf(MPI_COMM_WORLD,stderr,"\rProgress: %.2f%%",(double)nz/nnzwanted*100));
647 : }
648 12 : if (progress) PetscCall(PetscFPrintf(MPI_COMM_WORLD,stderr,"\r"));
649 :
650 12 : PetscCall(PetscFree2(iwork,swork));
651 12 : PetscCall(PetscRandomDestroy(&rand));
652 :
653 12 : PetscCall(XMatConvert(XA,A));
654 12 : PetscCall(XMatDestroy(&XA));
655 12 : PetscCall(MatGetInfo(A,MAT_GLOBAL_SUM,&Ainfo));
656 12 : 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)));
657 :
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
662 :
663 : /* scale by 0.5 the lower Q rows of A */
664 12 : PetscCall(MatCreateVecs(A,NULL,&vomega));
665 12 : PetscCall(VecSet(vomega,1.0));
666 652 : for (i=PetscMax(P,Istart);i<Iend;i++) {
667 640 : PetscCall(VecSetValue(vomega,i,0.5,INSERT_VALUES));
668 : }
669 12 : PetscCall(VecAssemblyBegin(vomega));
670 12 : PetscCall(VecAssemblyEnd(vomega));
671 12 : PetscCall(MatDiagonalScale(A,vomega,NULL));
672 :
673 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
674 : Create the signature
675 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
676 :
677 12 : PetscCall(VecSet(vomega,1.0));
678 12 : PetscCall(VecGetOwnershipRange(vomega,&Istart,&Iend));
679 1532 : for (i=Istart;i<Iend;i++) {
680 1520 : if (i>=P) PetscCall(VecSetValue(vomega,i,-1.0,INSERT_VALUES));
681 : }
682 12 : PetscCall(VecAssemblyBegin(vomega));
683 12 : PetscCall(VecAssemblyEnd(vomega));
684 :
685 12 : PetscCall(MatGetType(A,&Atype));
686 12 : PetscCall(MatCreate(PETSC_COMM_WORLD,&Omega));
687 12 : PetscCall(MatSetSizes(Omega,PETSC_DECIDE,PETSC_DECIDE,P+Q,P+Q));
688 12 : PetscCall(MatSetType(Omega,Atype));
689 12 : PetscCall(MatDiagonalSet(Omega,vomega,INSERT_VALUES));
690 :
691 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
692 : Create the singular value solver and set various options
693 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
694 :
695 12 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
696 :
697 12 : PetscCall(SVDSetOperators(svd,A,NULL));
698 12 : PetscCall(SVDSetSignature(svd,vomega));
699 12 : PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
700 :
701 12 : PetscCall(SVDSetFromOptions(svd));
702 :
703 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
704 : Solve the problem, display solution
705 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
706 :
707 12 : PetscCall(MatCreateVecs(A,&v,&u));
708 12 : PetscCall(SVDSolve(svd));
709 :
710 : /* show detailed info unless -terse option is given by user */
711 12 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
712 12 : if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
713 : else {
714 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
715 0 : PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
716 0 : PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
717 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
718 : }
719 :
720 : /* check orthogonality */
721 12 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-skiporth",&skiporth,NULL));
722 12 : PetscCall(SVDGetConverged(svd,&nconv));
723 12 : PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
724 12 : nconv = PetscMin(nconv,nsv);
725 12 : if (nconv>0 && !skiporth) {
726 12 : PetscCall(SVDGetTolerances(svd,&tol,NULL));
727 12 : PetscCall(VecDuplicateVecs(u,nconv,&U));
728 12 : PetscCall(VecDuplicateVecs(v,nconv,&V));
729 72 : for (i=0;i<nconv;i++) PetscCall(SVDGetSingularTriplet(svd,i,NULL,U[i],V[i]));
730 12 : PetscCall(VecCheckOrthonormality(U,nconv,NULL,nconv,Omega,NULL,&lev1));
731 12 : PetscCall(VecCheckOrthonormality(V,nconv,NULL,nconv,NULL,NULL,&lev2));
732 12 : if (lev1+lev2<20*tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n"));
733 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g (U) %g (V)\n",(double)lev1,(double)lev2));
734 12 : PetscCall(VecDestroyVecs(nconv,&U));
735 12 : PetscCall(VecDestroyVecs(nconv,&V));
736 : }
737 12 : PetscCall(VecDestroy(&u));
738 12 : PetscCall(VecDestroy(&v));
739 :
740 : /* free work space */
741 12 : PetscCall(SVDDestroy(&svd));
742 12 : PetscCall(MatDestroy(&Omega));
743 12 : PetscCall(MatDestroy(&A));
744 12 : PetscCall(VecDestroy(&vomega));
745 12 : PetscCall(SlepcFinalize());
746 : return 0;
747 : }
748 :
749 : /*TEST
750 :
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
765 :
766 : TEST*/
|
