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 : Private DS routines
12 : */
13 :
14 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
15 : #include <slepcblaslapack.h>
16 :
17 16867 : PetscErrorCode DSAllocateMat_Private(DS ds,DSMatType m)
18 : {
19 16867 : PetscInt n,d,rows=0,cols;
20 16867 : PetscBool ispep,isnep;
21 :
22 16867 : PetscFunctionBegin;
23 16867 : n = ds->ld;
24 16867 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,DSPEP,&ispep));
25 16867 : if (ispep) {
26 1136 : if (m==DS_MAT_A || m==DS_MAT_B || m==DS_MAT_W || m==DS_MAT_U || m==DS_MAT_X || m==DS_MAT_Y) {
27 1094 : PetscCall(DSPEPGetDegree(ds,&d));
28 1094 : n = d*ds->ld;
29 : }
30 : } else {
31 15731 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,DSNEP,&isnep));
32 15731 : if (isnep) {
33 635 : if (m==DS_MAT_Q || m==DS_MAT_Z || m==DS_MAT_U || m==DS_MAT_V || m==DS_MAT_X || m==DS_MAT_Y) {
34 20 : PetscCall(DSNEPGetMinimality(ds,&d));
35 20 : n = d*ds->ld;
36 : }
37 : }
38 : }
39 16867 : cols = n;
40 :
41 16867 : switch (m) {
42 : case DS_MAT_T:
43 : cols = PetscDefined(USE_COMPLEX)? 2: 3;
44 : rows = n;
45 : break;
46 154 : case DS_MAT_D:
47 154 : cols = 1;
48 154 : rows = n;
49 154 : break;
50 3561 : case DS_MAT_X:
51 3561 : rows = ds->ld;
52 3561 : break;
53 81 : case DS_MAT_Y:
54 81 : rows = ds->ld;
55 81 : break;
56 12497 : default:
57 12497 : rows = n;
58 : }
59 16867 : if (ds->omat[m]) PetscCall(MatZeroEntries(ds->omat[m]));
60 : else {
61 4201 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]));
62 : }
63 16867 : PetscFunctionReturn(PETSC_SUCCESS);
64 : }
65 :
66 41 : PetscErrorCode DSReallocateMat_Private(DS ds,DSMatType mat,PetscInt ld)
67 : {
68 41 : Mat M;
69 41 : PetscInt i,m,n,rows,cols;
70 41 : PetscScalar *dst;
71 41 : PetscReal *rdst;
72 41 : const PetscScalar *src;
73 41 : const PetscReal *rsrc;
74 :
75 41 : PetscFunctionBegin;
76 41 : PetscCall(MatGetSize(ds->omat[mat],&m,&n));
77 41 : rows = ld;
78 41 : cols = m==n? ld: n;
79 41 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&M));
80 41 : PetscCall(MatDenseGetArrayRead(ds->omat[mat],&src));
81 41 : PetscCall(MatDenseGetArray(M,&dst));
82 : if (mat==DS_MAT_T && PetscDefined(USE_COMPLEX)) {
83 : rsrc = (const PetscReal*)src;
84 : rdst = (PetscReal*)dst;
85 : for (i=0;i<3;i++) PetscCall(PetscArraycpy(rdst+i*ld,rsrc+i*ds->ld,ds->ld));
86 : } else {
87 387 : for (i=0;i<n;i++) PetscCall(PetscArraycpy(dst+i*ld,src+i*ds->ld,ds->ld));
88 : }
89 41 : PetscCall(MatDenseRestoreArrayRead(ds->omat[mat],&src));
90 41 : PetscCall(MatDenseRestoreArray(M,&dst));
91 41 : PetscCall(MatDestroy(&ds->omat[mat]));
92 41 : ds->omat[mat] = M;
93 41 : PetscFunctionReturn(PETSC_SUCCESS);
94 : }
95 :
96 89788 : PetscErrorCode DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
97 : {
98 89788 : PetscFunctionBegin;
99 89788 : if (s>ds->lwork) {
100 1549 : PetscCall(PetscFree(ds->work));
101 1549 : PetscCall(PetscMalloc1(s,&ds->work));
102 1549 : ds->lwork = s;
103 : }
104 89788 : if (r>ds->lrwork) {
105 579 : PetscCall(PetscFree(ds->rwork));
106 579 : PetscCall(PetscMalloc1(r,&ds->rwork));
107 579 : ds->lrwork = r;
108 : }
109 89788 : if (i>ds->liwork) {
110 621 : PetscCall(PetscFree(ds->iwork));
111 621 : PetscCall(PetscMalloc1(i,&ds->iwork));
112 621 : ds->liwork = i;
113 : }
114 89788 : PetscFunctionReturn(PETSC_SUCCESS);
115 : }
116 :
117 : /*@
118 : DSViewMat - Prints one of the internal DS matrices.
119 :
120 : Collective
121 :
122 : Input Parameters:
123 : + ds - the direct solver context
124 : . viewer - visualization context
125 : - m - matrix to display
126 :
127 : Note:
128 : Works only for ascii viewers. Set the viewer in Matlab format if
129 : want to paste into Matlab.
130 :
131 : Level: developer
132 :
133 : .seealso: DSView()
134 : @*/
135 30 : PetscErrorCode DSViewMat(DS ds,PetscViewer viewer,DSMatType m)
136 : {
137 30 : PetscInt i,j,rows,cols;
138 30 : const PetscScalar *M=NULL,*v;
139 30 : PetscViewerFormat format;
140 : #if defined(PETSC_USE_COMPLEX)
141 : PetscBool allreal = PETSC_TRUE;
142 : const PetscReal *vr;
143 : #endif
144 :
145 30 : PetscFunctionBegin;
146 30 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
147 90 : PetscValidLogicalCollectiveEnum(ds,m,3);
148 30 : DSCheckValidMat(ds,m,3);
149 30 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer));
150 30 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
151 30 : PetscCheckSameComm(ds,1,viewer,2);
152 30 : PetscCall(PetscViewerGetFormat(viewer,&format));
153 30 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
154 30 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
155 30 : PetscCall(DSMatGetSize(ds,m,&rows,&cols));
156 30 : PetscCall(MatDenseGetArrayRead(ds->omat[m],&M));
157 : #if defined(PETSC_USE_COMPLEX)
158 : /* determine if matrix has all real values */
159 : if (m!=DS_MAT_T && m!=DS_MAT_D) {
160 : /* determine if matrix has all real values */
161 : v = M;
162 : for (i=0;i<rows;i++)
163 : for (j=0;j<cols;j++)
164 : if (PetscImaginaryPart(v[i+j*ds->ld])) { allreal = PETSC_FALSE; break; }
165 : }
166 : if (m==DS_MAT_T) cols=3;
167 : #endif
168 30 : if (format == PETSC_VIEWER_ASCII_MATLAB) {
169 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
170 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s = [\n",DSMatName[m]));
171 30 : } else PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[m]));
172 :
173 420 : for (i=0;i<rows;i++) {
174 390 : v = M+i;
175 : #if defined(PETSC_USE_COMPLEX)
176 : vr = (const PetscReal*)M+i; /* handle compact storage, 2nd column is in imaginary part of 1st column */
177 : #endif
178 4185 : for (j=0;j<cols;j++) {
179 : #if defined(PETSC_USE_COMPLEX)
180 : if (m!=DS_MAT_T && m!=DS_MAT_D) {
181 : if (allreal) PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v)));
182 : else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v)));
183 : } else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*vr));
184 : #else
185 3795 : PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v));
186 : #endif
187 3795 : v += ds->ld;
188 : #if defined(PETSC_USE_COMPLEX)
189 : if (m==DS_MAT_T) vr += ds->ld;
190 : #endif
191 : }
192 390 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
193 : }
194 30 : PetscCall(MatDenseRestoreArrayRead(ds->omat[m],&M));
195 :
196 30 : if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
197 30 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
198 30 : PetscCall(PetscViewerFlush(viewer));
199 30 : PetscFunctionReturn(PETSC_SUCCESS);
200 : }
201 :
202 19539 : PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
203 : {
204 19539 : PetscScalar re,im,wi0;
205 19539 : PetscInt n,i,j,result,tmp1,tmp2=0,d=1;
206 :
207 19539 : PetscFunctionBegin;
208 19539 : n = ds->t; /* sort only first t pairs if truncated */
209 : /* insertion sort */
210 19539 : i=ds->l+1;
211 : #if !defined(PETSC_USE_COMPLEX)
212 19539 : if (wi && wi[perm[i-1]]!=0.0) i++; /* initial value is complex */
213 : #else
214 : if (isghiep && PetscImaginaryPart(wr[perm[i-1]])!=0.0) i++;
215 : #endif
216 267773 : for (;i<n;i+=d) {
217 248234 : re = wr[perm[i]];
218 248234 : if (wi) im = wi[perm[i]];
219 : else im = 0.0;
220 248234 : tmp1 = perm[i];
221 : #if !defined(PETSC_USE_COMPLEX)
222 248234 : if (im!=0.0) { d = 2; tmp2 = perm[i+1]; }
223 : else d = 1;
224 : #else
225 : if (isghiep && PetscImaginaryPart(re)!=0.0) { d = 2; tmp2 = perm[i+1]; }
226 : else d = 1;
227 : #endif
228 248234 : j = i-1;
229 248234 : if (wi) wi0 = wi[perm[j]];
230 : else wi0 = 0.0;
231 248234 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
232 1145186 : while (result<0 && j>=ds->l) {
233 1067804 : perm[j+d] = perm[j];
234 1067804 : j--;
235 : #if !defined(PETSC_USE_COMPLEX)
236 1067804 : if (wi && wi[perm[j+1]]!=0)
237 : #else
238 : if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
239 : #endif
240 39323 : { perm[j+d] = perm[j]; j--; }
241 :
242 1067804 : if (j>=ds->l) {
243 990422 : if (wi) wi0 = wi[perm[j]];
244 : else wi0 = 0.0;
245 2306460 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
246 : }
247 : }
248 248234 : perm[j+1] = tmp1;
249 248234 : if (d==2) perm[j+2] = tmp2;
250 : }
251 19539 : PetscFunctionReturn(PETSC_SUCCESS);
252 : }
253 :
254 7002 : PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
255 : {
256 7002 : PetscScalar re;
257 7002 : PetscInt i,j,result,tmp,l,n;
258 :
259 7002 : PetscFunctionBegin;
260 7002 : n = ds->t; /* sort only first t pairs if truncated */
261 7002 : l = ds->l;
262 : /* insertion sort */
263 100480 : for (i=l+1;i<n;i++) {
264 93478 : re = eig[perm[i]];
265 93478 : j = i-1;
266 93478 : PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
267 526663 : while (result<0 && j>=l) {
268 497850 : tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
269 1089178 : if (j>=l) PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
270 : }
271 : }
272 7002 : PetscFunctionReturn(PETSC_SUCCESS);
273 : }
274 :
275 : /*
276 : Permute comumns [istart..iend-1] of [mat] according to perm. Columns have length n
277 : */
278 22987 : PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat,PetscInt *perm)
279 : {
280 22987 : PetscInt i,j,k,p,ld;
281 22987 : PetscScalar *Q,rtmp;
282 :
283 22987 : PetscFunctionBegin;
284 22987 : ld = ds->ld;
285 22987 : PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
286 346833 : for (i=istart;i<iend;i++) {
287 323846 : p = perm[i];
288 323846 : if (p != i) {
289 144179 : j = i + 1;
290 1202040 : while (perm[j] != i) j++;
291 144179 : perm[j] = p; perm[i] = i;
292 : /* swap columns i and j */
293 4362175 : for (k=0;k<n;k++) {
294 4217996 : rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
295 : }
296 : }
297 : }
298 22987 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
299 22987 : PetscFunctionReturn(PETSC_SUCCESS);
300 : }
301 :
302 : /*
303 : The same as DSPermuteColumns_Private but for two matrices [mat1] and [mat2]
304 : */
305 267 : PetscErrorCode DSPermuteColumnsTwo_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
306 : {
307 267 : PetscInt i,j,k,p,ld;
308 267 : PetscScalar *Q,*Z,rtmp,rtmp2;
309 :
310 267 : PetscFunctionBegin;
311 267 : ld = ds->ld;
312 267 : PetscCall(MatDenseGetArray(ds->omat[mat1],&Q));
313 267 : PetscCall(MatDenseGetArray(ds->omat[mat2],&Z));
314 5813 : for (i=istart;i<iend;i++) {
315 5546 : p = perm[i];
316 5546 : if (p != i) {
317 4021 : j = i + 1;
318 38032 : while (perm[j] != i) j++;
319 4021 : perm[j] = p; perm[i] = i;
320 : /* swap columns i and j */
321 56518 : for (k=0;k<n;k++) {
322 52497 : rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
323 52497 : rtmp2 = Z[k+p*ld]; Z[k+p*ld] = Z[k+i*ld]; Z[k+i*ld] = rtmp2;
324 : }
325 : }
326 : }
327 267 : PetscCall(MatDenseRestoreArray(ds->omat[mat1],&Q));
328 267 : PetscCall(MatDenseRestoreArray(ds->omat[mat2],&Z));
329 267 : PetscFunctionReturn(PETSC_SUCCESS);
330 : }
331 :
332 : /*
333 : Permute rows [istart..iend-1] of [mat] according to perm. Rows have length m
334 : */
335 0 : PetscErrorCode DSPermuteRows_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt m,DSMatType mat,PetscInt *perm)
336 : {
337 0 : PetscInt i,j,k,p,ld;
338 0 : PetscScalar *Q,rtmp;
339 :
340 0 : PetscFunctionBegin;
341 0 : ld = ds->ld;
342 0 : PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
343 0 : for (i=istart;i<iend;i++) {
344 0 : p = perm[i];
345 0 : if (p != i) {
346 0 : j = i + 1;
347 0 : while (perm[j] != i) j++;
348 0 : perm[j] = p; perm[i] = i;
349 : /* swap rows i and j */
350 0 : for (k=0;k<m;k++) {
351 0 : rtmp = Q[p+k*ld]; Q[p+k*ld] = Q[i+k*ld]; Q[i+k*ld] = rtmp;
352 : }
353 : }
354 : }
355 0 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
356 0 : PetscFunctionReturn(PETSC_SUCCESS);
357 : }
358 :
359 : /*
360 : Permute columns [istart..iend-1] of [mat1] and [mat2] according to perm.
361 : Columns of [mat1] have length n, columns of [mat2] have length m
362 : */
363 2438 : PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,PetscInt m,DSMatType mat1,DSMatType mat2,PetscInt *perm)
364 : {
365 2438 : PetscInt i,j,k,p,ld;
366 2438 : PetscScalar *U,*V,rtmp;
367 :
368 2438 : PetscFunctionBegin;
369 2438 : ld = ds->ld;
370 2438 : PetscCall(MatDenseGetArray(ds->omat[mat1],&U));
371 2438 : PetscCall(MatDenseGetArray(ds->omat[mat2],&V));
372 32305 : for (i=istart;i<iend;i++) {
373 29867 : p = perm[i];
374 29867 : if (p != i) {
375 6428 : j = i + 1;
376 81491 : while (perm[j] != i) j++;
377 6428 : perm[j] = p; perm[i] = i;
378 : /* swap columns i and j of U */
379 498951 : for (k=0;k<n;k++) {
380 492523 : rtmp = U[k+p*ld]; U[k+p*ld] = U[k+i*ld]; U[k+i*ld] = rtmp;
381 : }
382 : /* swap columns i and j of V */
383 297132 : for (k=0;k<m;k++) {
384 290704 : rtmp = V[k+p*ld]; V[k+p*ld] = V[k+i*ld]; V[k+i*ld] = rtmp;
385 : }
386 : }
387 : }
388 2438 : PetscCall(MatDenseRestoreArray(ds->omat[mat1],&U));
389 2438 : PetscCall(MatDenseRestoreArray(ds->omat[mat2],&V));
390 2438 : PetscFunctionReturn(PETSC_SUCCESS);
391 : }
392 :
393 : /*@
394 : DSSetIdentity - Copy the identity (a diagonal matrix with ones) on the
395 : active part of a matrix.
396 :
397 : Logically Collective
398 :
399 : Input Parameters:
400 : + ds - the direct solver context
401 : - mat - the matrix to modify
402 :
403 : Level: intermediate
404 :
405 : .seealso: DSGetMat()
406 : @*/
407 10267 : PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
408 : {
409 10267 : PetscScalar *A;
410 10267 : PetscInt i,ld,n,l;
411 :
412 10267 : PetscFunctionBegin;
413 10267 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
414 30801 : PetscValidLogicalCollectiveEnum(ds,mat,2);
415 10267 : DSCheckValidMat(ds,mat,2);
416 :
417 10267 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
418 10267 : PetscCall(DSGetLeadingDimension(ds,&ld));
419 10267 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
420 10267 : PetscCall(MatDenseGetArray(ds->omat[mat],&A));
421 10267 : PetscCall(PetscArrayzero(A+l*ld,ld*(n-l)));
422 144351 : for (i=l;i<n;i++) A[i+i*ld] = 1.0;
423 10267 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
424 10267 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
425 10267 : PetscFunctionReturn(PETSC_SUCCESS);
426 : }
427 :
428 : /*@
429 : DSOrthogonalize - Orthogonalize the columns of a matrix.
430 :
431 : Logically Collective
432 :
433 : Input Parameters:
434 : + ds - the direct solver context
435 : . mat - a matrix
436 : - cols - number of columns to orthogonalize (starting from column zero)
437 :
438 : Output Parameter:
439 : . lindcols - (optional) number of linearly independent columns of the matrix
440 :
441 : Level: developer
442 :
443 : .seealso: DSPseudoOrthogonalize()
444 : @*/
445 1614 : PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
446 : {
447 1614 : PetscInt n,l,ld;
448 1614 : PetscBLASInt ld_,rA,cA,info,ltau,lw;
449 1614 : PetscScalar *A,*tau,*w,saux,dummy;
450 :
451 1614 : PetscFunctionBegin;
452 1614 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
453 1614 : DSCheckAlloc(ds,1);
454 4842 : PetscValidLogicalCollectiveEnum(ds,mat,2);
455 1614 : DSCheckValidMat(ds,mat,2);
456 4842 : PetscValidLogicalCollectiveInt(ds,cols,3);
457 :
458 1614 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
459 1614 : PetscCall(DSGetLeadingDimension(ds,&ld));
460 1614 : n = n - l;
461 1614 : PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
462 1614 : if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
463 :
464 1614 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
465 1614 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
466 1614 : PetscCall(MatDenseGetArray(ds->omat[mat],&A));
467 1614 : PetscCall(PetscBLASIntCast(PetscMin(cols,n),<au));
468 1614 : PetscCall(PetscBLASIntCast(ld,&ld_));
469 1614 : PetscCall(PetscBLASIntCast(n,&rA));
470 1614 : PetscCall(PetscBLASIntCast(cols,&cA));
471 1614 : lw = -1;
472 1614 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
473 1614 : SlepcCheckLapackInfo("geqrf",info);
474 1614 : lw = (PetscBLASInt)PetscRealPart(saux);
475 1614 : PetscCall(DSAllocateWork_Private(ds,lw+ltau,0,0));
476 1614 : tau = ds->work;
477 1614 : w = &tau[ltau];
478 1614 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
479 1614 : SlepcCheckLapackInfo("geqrf",info);
480 1614 : PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,<au,<au,&A[ld*l+l],&ld_,tau,w,&lw,&info));
481 1614 : SlepcCheckLapackInfo("orgqr",info);
482 1614 : if (lindcols) *lindcols = ltau;
483 1614 : PetscCall(PetscFPTrapPop());
484 1614 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
485 1614 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
486 1614 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
487 1614 : PetscFunctionReturn(PETSC_SUCCESS);
488 : }
489 :
490 : /*
491 : Compute C <- a*A*B + b*C, where
492 : ldC, the leading dimension of C,
493 : ldA, the leading dimension of A,
494 : rA, cA, rows and columns of A,
495 : At, if true use the transpose of A instead,
496 : ldB, the leading dimension of B,
497 : rB, cB, rows and columns of B,
498 : Bt, if true use the transpose of B instead
499 : */
500 55 : static PetscErrorCode SlepcMatDenseMult(PetscScalar *C,PetscInt _ldC,PetscScalar b,PetscScalar a,const PetscScalar *A,PetscInt _ldA,PetscInt rA,PetscInt cA,PetscBool At,const PetscScalar *B,PetscInt _ldB,PetscInt rB,PetscInt cB,PetscBool Bt)
501 : {
502 55 : PetscInt tmp;
503 55 : PetscBLASInt m, n, k, ldA, ldB, ldC;
504 55 : const char *N = "N", *T = "C", *qA = N, *qB = N;
505 :
506 55 : PetscFunctionBegin;
507 55 : if ((rA == 0) || (cB == 0)) PetscFunctionReturn(PETSC_SUCCESS);
508 55 : PetscAssertPointer(C,1);
509 55 : PetscAssertPointer(A,5);
510 55 : PetscAssertPointer(B,10);
511 55 : PetscCall(PetscBLASIntCast(_ldA,&ldA));
512 55 : PetscCall(PetscBLASIntCast(_ldB,&ldB));
513 55 : PetscCall(PetscBLASIntCast(_ldC,&ldC));
514 :
515 : /* Transpose if needed */
516 55 : if (At) tmp = rA, rA = cA, cA = tmp, qA = T;
517 55 : if (Bt) tmp = rB, rB = cB, cB = tmp, qB = T;
518 :
519 : /* Check size */
520 55 : PetscCheck(cA==rB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix dimensions do not match");
521 :
522 : /* Do stub */
523 55 : if ((rA == 1) && (cA == 1) && (cB == 1)) {
524 0 : if (!At && !Bt) *C = *A * *B;
525 0 : else if (At && !Bt) *C = PetscConj(*A) * *B;
526 0 : else if (!At && Bt) *C = *A * PetscConj(*B);
527 0 : else *C = PetscConj(*A) * PetscConj(*B);
528 0 : m = n = k = 1;
529 : } else {
530 55 : PetscCall(PetscBLASIntCast(rA,&m));
531 55 : PetscCall(PetscBLASIntCast(cB,&n));
532 55 : PetscCall(PetscBLASIntCast(cA,&k));
533 55 : PetscCallBLAS("BLASgemm",BLASgemm_(qA,qB,&m,&n,&k,&a,(PetscScalar*)A,&ldA,(PetscScalar*)B,&ldB,&b,C,&ldC));
534 : }
535 :
536 55 : PetscCall(PetscLogFlops(2.0*m*n*k));
537 55 : PetscFunctionReturn(PETSC_SUCCESS);
538 : }
539 :
540 : /*@
541 : DSPseudoOrthogonalize - Orthogonalize the columns of a matrix with Modified
542 : Gram-Schmidt in an indefinite inner product space defined by a signature.
543 :
544 : Logically Collective
545 :
546 : Input Parameters:
547 : + ds - the direct solver context
548 : . mat - the matrix
549 : . cols - number of columns to orthogonalize (starting from column zero)
550 : - s - the signature that defines the inner product
551 :
552 : Output Parameters:
553 : + lindcols - (optional) linearly independent columns of the matrix
554 : - ns - (optional) the new signature of the vectors
555 :
556 : Note:
557 : After the call the matrix satisfies A'*s*A = ns.
558 :
559 : Level: developer
560 :
561 : .seealso: DSOrthogonalize()
562 : @*/
563 1 : PetscErrorCode DSPseudoOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscReal s[],PetscInt *lindcols,PetscReal ns[])
564 : {
565 1 : PetscInt i,j,k,l,n,ld;
566 1 : PetscBLASInt info,one=1,zero=0,rA_,ld_;
567 1 : PetscScalar *A,*A_,*m,*h,nr0;
568 1 : PetscReal nr_o,nr,nr_abs,*ns_,done=1.0;
569 :
570 1 : PetscFunctionBegin;
571 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
572 1 : DSCheckAlloc(ds,1);
573 3 : PetscValidLogicalCollectiveEnum(ds,mat,2);
574 1 : DSCheckValidMat(ds,mat,2);
575 3 : PetscValidLogicalCollectiveInt(ds,cols,3);
576 1 : PetscAssertPointer(s,4);
577 1 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
578 1 : PetscCall(DSGetLeadingDimension(ds,&ld));
579 1 : n = n - l;
580 1 : PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
581 1 : if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
582 1 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
583 1 : PetscCall(PetscBLASIntCast(n,&rA_));
584 1 : PetscCall(PetscBLASIntCast(ld,&ld_));
585 1 : PetscCall(MatDenseGetArray(ds->omat[mat],&A_));
586 1 : A = &A_[ld*l+l];
587 2 : PetscCall(DSAllocateWork_Private(ds,n+cols,ns?0:cols,0));
588 1 : m = ds->work;
589 1 : h = &m[n];
590 1 : ns_ = ns ? ns : ds->rwork;
591 11 : for (i=0; i<cols; i++) {
592 : /* m <- diag(s)*A[i] */
593 110 : for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
594 : /* nr_o <- mynorm(A[i]'*m), mynorm(x) = sign(x)*sqrt(|x|) */
595 10 : PetscCall(SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
596 10 : nr = nr_o = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
597 16 : for (j=0; j<3 && i>0; j++) {
598 : /* h <- A[0:i-1]'*m */
599 15 : PetscCall(SlepcMatDenseMult(h,i,0.0,1.0,A,ld,n,i,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
600 : /* h <- diag(ns)*h */
601 99 : for (k=0; k<i; k++) h[k] *= ns_[k];
602 : /* A[i] <- A[i] - A[0:i-1]*h */
603 15 : PetscCall(SlepcMatDenseMult(&A[ld*i],ld,1.0,-1.0,A,ld,n,i,PETSC_FALSE,h,i,i,1,PETSC_FALSE));
604 : /* m <- diag(s)*A[i] */
605 165 : for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
606 : /* nr_o <- mynorm(A[i]'*m) */
607 15 : PetscCall(SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
608 15 : nr = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
609 15 : PetscCheck(PetscAbs(nr)>PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linear dependency detected");
610 15 : if (PetscAbs(nr) > 0.7*PetscAbs(nr_o)) break;
611 6 : nr_o = nr;
612 : }
613 10 : ns_[i] = PetscSign(nr);
614 : /* A[i] <- A[i]/|nr| */
615 10 : nr_abs = PetscAbs(nr);
616 10 : PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nr_abs,&done,&rA_,&one,A+i*ld,&ld_,&info));
617 10 : SlepcCheckLapackInfo("lascl",info);
618 : }
619 1 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A_));
620 1 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
621 1 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
622 1 : if (lindcols) *lindcols = cols;
623 1 : PetscFunctionReturn(PETSC_SUCCESS);
624 : }
|
