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 7700 : PetscErrorCode DSAllocateMat_Private(DS ds,DSMatType m)
18 : {
19 7700 : PetscInt n,d,rows=0,cols;
20 7700 : PetscBool ispep,isnep;
21 :
22 7700 : PetscFunctionBegin;
23 7700 : n = ds->ld;
24 7700 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,DSPEP,&ispep));
25 7700 : if (ispep) {
26 1389 : 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 1304 : PetscCall(DSPEPGetDegree(ds,&d));
28 1304 : n = d*ds->ld;
29 : }
30 : } else {
31 6311 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,DSNEP,&isnep));
32 6311 : if (isnep) {
33 775 : 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 152 : PetscCall(DSNEPGetMinimality(ds,&d));
35 152 : n = d*ds->ld;
36 : }
37 : }
38 : }
39 7700 : cols = n;
40 :
41 7700 : switch (m) {
42 : case DS_MAT_T:
43 : cols = PetscDefined(USE_COMPLEX)? 2: 3;
44 : rows = n;
45 : break;
46 140 : case DS_MAT_D:
47 140 : cols = 1;
48 140 : rows = n;
49 140 : break;
50 1216 : case DS_MAT_X:
51 1216 : rows = ds->ld;
52 1216 : break;
53 94 : case DS_MAT_Y:
54 94 : rows = ds->ld;
55 94 : break;
56 5708 : default:
57 5708 : rows = n;
58 : }
59 7700 : if (ds->omat[m]) PetscCall(MatZeroEntries(ds->omat[m]));
60 : else {
61 4467 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]));
62 : }
63 7700 : PetscFunctionReturn(PETSC_SUCCESS);
64 : }
65 :
66 30 : PetscErrorCode DSReallocateMat_Private(DS ds,DSMatType mat,PetscInt ld)
67 : {
68 30 : Mat M;
69 30 : PetscInt i,m,n,rows,cols;
70 30 : PetscScalar *dst;
71 30 : PetscReal *rdst;
72 30 : const PetscScalar *src;
73 30 : const PetscReal *rsrc;
74 :
75 30 : PetscFunctionBegin;
76 30 : PetscCall(MatGetSize(ds->omat[mat],&m,&n));
77 30 : rows = ld;
78 30 : cols = m==n? ld: n;
79 30 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&M));
80 30 : PetscCall(MatDenseGetArrayRead(ds->omat[mat],&src));
81 30 : PetscCall(MatDenseGetArray(M,&dst));
82 30 : if (mat==DS_MAT_T && PetscDefined(USE_COMPLEX)) {
83 11 : rsrc = (const PetscReal*)src;
84 11 : rdst = (PetscReal*)dst;
85 44 : for (i=0;i<3;i++) PetscCall(PetscArraycpy(rdst+i*ld,rsrc+i*ds->ld,ds->ld));
86 : } else {
87 242 : for (i=0;i<n;i++) PetscCall(PetscArraycpy(dst+i*ld,src+i*ds->ld,ds->ld));
88 : }
89 30 : PetscCall(MatDenseRestoreArrayRead(ds->omat[mat],&src));
90 30 : PetscCall(MatDenseRestoreArray(M,&dst));
91 30 : PetscCall(MatDestroy(&ds->omat[mat]));
92 30 : ds->omat[mat] = M;
93 30 : PetscFunctionReturn(PETSC_SUCCESS);
94 : }
95 :
96 79407 : PetscErrorCode DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
97 : {
98 79407 : PetscFunctionBegin;
99 79407 : if (s>ds->lwork) {
100 1611 : PetscCall(PetscFree(ds->work));
101 1611 : PetscCall(PetscMalloc1(s,&ds->work));
102 1611 : ds->lwork = s;
103 : }
104 79407 : if (r>ds->lrwork) {
105 1354 : PetscCall(PetscFree(ds->rwork));
106 1354 : PetscCall(PetscMalloc1(r,&ds->rwork));
107 1354 : ds->lrwork = r;
108 : }
109 79407 : if (i>ds->liwork) {
110 651 : PetscCall(PetscFree(ds->iwork));
111 651 : PetscCall(PetscMalloc1(i,&ds->iwork));
112 651 : ds->liwork = i;
113 : }
114 79407 : 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 30 : PetscBool allreal = PETSC_TRUE;
142 30 : 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 30 : if (m!=DS_MAT_T && m!=DS_MAT_D) {
160 : /* determine if matrix has all real values */
161 29 : v = M;
162 404 : for (i=0;i<rows;i++)
163 4155 : for (j=0;j<cols;j++)
164 3780 : if (PetscImaginaryPart(v[i+j*ds->ld])) { allreal = PETSC_FALSE; break; }
165 : }
166 30 : 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 390 : 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 3795 : if (m!=DS_MAT_T && m!=DS_MAT_D) {
181 3780 : if (allreal) PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v)));
182 0 : else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v)));
183 15 : } else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*vr));
184 : #else
185 : PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v));
186 : #endif
187 3795 : v += ds->ld;
188 : #if defined(PETSC_USE_COMPLEX)
189 3795 : 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 19032 : PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
203 : {
204 19032 : PetscScalar re,im,wi0;
205 19032 : PetscInt n,i,j,result,tmp1,tmp2=0,d=1;
206 :
207 19032 : PetscFunctionBegin;
208 19032 : n = ds->t; /* sort only first t pairs if truncated */
209 : /* insertion sort */
210 19032 : i=ds->l+1;
211 : #if !defined(PETSC_USE_COMPLEX)
212 : if (wi && wi[perm[i-1]]!=0.0) i++; /* initial value is complex */
213 : #else
214 19032 : if (isghiep && PetscImaginaryPart(wr[perm[i-1]])!=0.0) i++;
215 : #endif
216 221965 : for (;i<n;i+=d) {
217 202933 : re = wr[perm[i]];
218 202933 : if (wi) im = wi[perm[i]];
219 : else im = 0.0;
220 202933 : tmp1 = perm[i];
221 : #if !defined(PETSC_USE_COMPLEX)
222 : if (im!=0.0) { d = 2; tmp2 = perm[i+1]; }
223 : else d = 1;
224 : #else
225 202933 : if (isghiep && PetscImaginaryPart(re)!=0.0) { d = 2; tmp2 = perm[i+1]; }
226 : else d = 1;
227 : #endif
228 202933 : j = i-1;
229 202933 : if (wi) wi0 = wi[perm[j]];
230 : else wi0 = 0.0;
231 202933 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
232 622030 : while (result<0 && j>=ds->l) {
233 568800 : perm[j+d] = perm[j];
234 568800 : j--;
235 : #if !defined(PETSC_USE_COMPLEX)
236 : if (wi && wi[perm[j+1]]!=0)
237 : #else
238 568800 : if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
239 : #endif
240 3911 : { perm[j+d] = perm[j]; j--; }
241 :
242 568800 : if (j>=ds->l) {
243 515570 : if (wi) wi0 = wi[perm[j]];
244 : else wi0 = 0.0;
245 1287303 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
246 : }
247 : }
248 202933 : perm[j+1] = tmp1;
249 202933 : if (d==2) perm[j+2] = tmp2;
250 : }
251 19032 : PetscFunctionReturn(PETSC_SUCCESS);
252 : }
253 :
254 5492 : PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
255 : {
256 5492 : PetscScalar re;
257 5492 : PetscInt i,j,result,tmp,l,n;
258 :
259 5492 : PetscFunctionBegin;
260 5492 : n = ds->t; /* sort only first t pairs if truncated */
261 5492 : l = ds->l;
262 : /* insertion sort */
263 73768 : for (i=l+1;i<n;i++) {
264 68276 : re = eig[perm[i]];
265 68276 : j = i-1;
266 68276 : PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
267 525598 : while (result<0 && j>=l) {
268 498369 : tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
269 1065014 : if (j>=l) PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
270 : }
271 : }
272 5492 : PetscFunctionReturn(PETSC_SUCCESS);
273 : }
274 :
275 : /*
276 : Permute comumns [istart..iend-1] of [mat] according to perm. Columns have length n
277 : */
278 20369 : PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat,PetscInt *perm)
279 : {
280 20369 : PetscInt i,j,k,p,ld;
281 20369 : PetscScalar *Q,rtmp;
282 :
283 20369 : PetscFunctionBegin;
284 20369 : ld = ds->ld;
285 20369 : PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
286 263344 : for (i=istart;i<iend;i++) {
287 242975 : p = perm[i];
288 242975 : if (p != i) {
289 72605 : j = i + 1;
290 571246 : while (perm[j] != i) j++;
291 72605 : perm[j] = p; perm[i] = i;
292 : /* swap columns i and j */
293 1982440 : for (k=0;k<n;k++) {
294 1909835 : rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
295 : }
296 : }
297 : }
298 20369 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
299 20369 : PetscFunctionReturn(PETSC_SUCCESS);
300 : }
301 :
302 : /*
303 : The same as DSPermuteColumns_Private but for two matrices [mat1] and [mat2]
304 : */
305 314 : PetscErrorCode DSPermuteColumnsTwo_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
306 : {
307 314 : PetscInt i,j,k,p,ld;
308 314 : PetscScalar *Q,*Z,rtmp,rtmp2;
309 :
310 314 : PetscFunctionBegin;
311 314 : ld = ds->ld;
312 314 : PetscCall(MatDenseGetArray(ds->omat[mat1],&Q));
313 314 : PetscCall(MatDenseGetArray(ds->omat[mat2],&Z));
314 7963 : for (i=istart;i<iend;i++) {
315 7649 : p = perm[i];
316 7649 : if (p != i) {
317 5395 : j = i + 1;
318 51762 : while (perm[j] != i) j++;
319 5395 : perm[j] = p; perm[i] = i;
320 : /* swap columns i and j */
321 95509 : for (k=0;k<n;k++) {
322 90114 : rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
323 90114 : rtmp2 = Z[k+p*ld]; Z[k+p*ld] = Z[k+i*ld]; Z[k+i*ld] = rtmp2;
324 : }
325 : }
326 : }
327 314 : PetscCall(MatDenseRestoreArray(ds->omat[mat1],&Q));
328 314 : PetscCall(MatDenseRestoreArray(ds->omat[mat2],&Z));
329 314 : 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 1885 : PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,PetscInt m,DSMatType mat1,DSMatType mat2,PetscInt *perm)
364 : {
365 1885 : PetscInt i,j,k,p,ld;
366 1885 : PetscScalar *U,*V,rtmp;
367 :
368 1885 : PetscFunctionBegin;
369 1885 : ld = ds->ld;
370 1885 : PetscCall(MatDenseGetArray(ds->omat[mat1],&U));
371 1885 : PetscCall(MatDenseGetArray(ds->omat[mat2],&V));
372 21428 : for (i=istart;i<iend;i++) {
373 19543 : p = perm[i];
374 19543 : if (p != i) {
375 4088 : j = i + 1;
376 66926 : while (perm[j] != i) j++;
377 4088 : perm[j] = p; perm[i] = i;
378 : /* swap columns i and j of U */
379 433052 : for (k=0;k<n;k++) {
380 428964 : 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 230520 : for (k=0;k<m;k++) {
384 226432 : rtmp = V[k+p*ld]; V[k+p*ld] = V[k+i*ld]; V[k+i*ld] = rtmp;
385 : }
386 : }
387 : }
388 1885 : PetscCall(MatDenseRestoreArray(ds->omat[mat1],&U));
389 1885 : PetscCall(MatDenseRestoreArray(ds->omat[mat2],&V));
390 1885 : 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 10110 : PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
408 : {
409 10110 : PetscScalar *A;
410 10110 : PetscInt i,ld,n,l;
411 :
412 10110 : PetscFunctionBegin;
413 10110 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
414 30330 : PetscValidLogicalCollectiveEnum(ds,mat,2);
415 10110 : DSCheckValidMat(ds,mat,2);
416 :
417 10110 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
418 10110 : PetscCall(DSGetLeadingDimension(ds,&ld));
419 10110 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
420 10110 : PetscCall(MatDenseGetArray(ds->omat[mat],&A));
421 10110 : PetscCall(PetscArrayzero(A+l*ld,ld*(n-l)));
422 137644 : for (i=l;i<n;i++) A[i+i*ld] = 1.0;
423 10110 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
424 10110 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
425 10110 : 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 1778 : PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
446 : {
447 1778 : PetscInt n,l,ld;
448 1778 : PetscBLASInt ld_,rA,cA,info,ltau,lw;
449 1778 : PetscScalar *A,*tau,*w,saux,dummy;
450 :
451 1778 : PetscFunctionBegin;
452 1778 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
453 1778 : DSCheckAlloc(ds,1);
454 5334 : PetscValidLogicalCollectiveEnum(ds,mat,2);
455 1778 : DSCheckValidMat(ds,mat,2);
456 5334 : PetscValidLogicalCollectiveInt(ds,cols,3);
457 :
458 1778 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
459 1778 : PetscCall(DSGetLeadingDimension(ds,&ld));
460 1778 : n = n - l;
461 1778 : PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
462 1778 : if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
463 :
464 1778 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
465 1778 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
466 1778 : PetscCall(MatDenseGetArray(ds->omat[mat],&A));
467 1778 : PetscCall(PetscBLASIntCast(PetscMin(cols,n),<au));
468 1778 : PetscCall(PetscBLASIntCast(ld,&ld_));
469 1778 : PetscCall(PetscBLASIntCast(n,&rA));
470 1778 : PetscCall(PetscBLASIntCast(cols,&cA));
471 1778 : lw = -1;
472 1778 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
473 1778 : SlepcCheckLapackInfo("geqrf",info);
474 1778 : lw = (PetscBLASInt)PetscRealPart(saux);
475 1778 : PetscCall(DSAllocateWork_Private(ds,lw+ltau,0,0));
476 1778 : tau = ds->work;
477 1778 : w = &tau[ltau];
478 1778 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
479 1778 : SlepcCheckLapackInfo("geqrf",info);
480 1778 : PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,<au,<au,&A[ld*l+l],&ld_,tau,w,&lw,&info));
481 1778 : SlepcCheckLapackInfo("orgqr",info);
482 1778 : if (lindcols) *lindcols = ltau;
483 1778 : PetscCall(PetscFPTrapPop());
484 1778 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
485 1778 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
486 1778 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
487 1778 : 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 : }
|
