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