Actual source code: dsghep.c
slepc-3.18.3 2023-03-24
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/dsimpl.h>
12: #include <slepcblaslapack.h>
14: PetscErrorCode DSAllocate_GHEP(DS ds,PetscInt ld)
15: {
16: DSAllocateMat_Private(ds,DS_MAT_A);
17: DSAllocateMat_Private(ds,DS_MAT_B);
18: DSAllocateMat_Private(ds,DS_MAT_Q);
19: PetscFree(ds->perm);
20: PetscMalloc1(ld,&ds->perm);
21: return 0;
22: }
24: PetscErrorCode DSView_GHEP(DS ds,PetscViewer viewer)
25: {
26: PetscViewerFormat format;
28: PetscViewerGetFormat(viewer,&format);
29: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
30: DSViewMat(ds,viewer,DS_MAT_A);
31: DSViewMat(ds,viewer,DS_MAT_B);
32: if (ds->state>DS_STATE_INTERMEDIATE) DSViewMat(ds,viewer,DS_MAT_Q);
33: if (ds->omat[DS_MAT_X]) DSViewMat(ds,viewer,DS_MAT_X);
34: return 0;
35: }
37: PetscErrorCode DSVectors_GHEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
38: {
39: PetscScalar *Z;
40: const PetscScalar *Q;
41: PetscInt ld = ds->ld;
44: switch (mat) {
45: case DS_MAT_X:
46: case DS_MAT_Y:
47: if (j) {
48: MatDenseGetArray(ds->omat[mat],&Z);
49: if (ds->state>=DS_STATE_CONDENSED) {
50: MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q);
51: PetscArraycpy(Z+(*j)*ld,Q+(*j)*ld,ld);
52: MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q);
53: } else {
54: PetscArrayzero(Z+(*j)*ld,ld);
55: Z[(*j)+(*j)*ld] = 1.0;
56: }
57: MatDenseRestoreArray(ds->omat[mat],&Z);
58: } else {
59: if (ds->state>=DS_STATE_CONDENSED) MatCopy(ds->omat[DS_MAT_Q],ds->omat[mat],SAME_NONZERO_PATTERN);
60: else DSSetIdentity(ds,mat);
61: }
62: break;
63: case DS_MAT_U:
64: case DS_MAT_V:
65: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
66: default:
67: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
68: }
69: return 0;
70: }
72: PetscErrorCode DSSort_GHEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
73: {
74: PetscInt n,l,i,*perm,ld=ds->ld;
75: PetscScalar *A;
77: if (!ds->sc) return 0;
78: n = ds->n;
79: l = ds->l;
80: MatDenseGetArray(ds->omat[DS_MAT_A],&A);
81: perm = ds->perm;
82: for (i=l;i<n;i++) wr[i] = A[i+i*ld];
83: if (rr) DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
84: else DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE);
85: for (i=l;i<n;i++) A[i+i*ld] = wr[perm[i]];
86: for (i=l;i<n;i++) wr[i] = A[i+i*ld];
87: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
88: DSPermuteColumns_Private(ds,l,n,n,DS_MAT_Q,perm);
89: return 0;
90: }
92: PetscErrorCode DSSolve_GHEP(DS ds,PetscScalar *wr,PetscScalar *wi)
93: {
94: PetscScalar *work,*A,*B,*Q;
95: PetscBLASInt itype = 1,*iwork,info,n1,liwork,ld,lrwork=0,lwork;
96: PetscInt off,i;
97: #if defined(PETSC_USE_COMPLEX)
98: PetscReal *rwork,*rr;
99: #endif
101: PetscBLASIntCast(ds->n-ds->l,&n1);
102: PetscBLASIntCast(ds->ld,&ld);
103: PetscBLASIntCast(5*ds->n+3,&liwork);
104: #if defined(PETSC_USE_COMPLEX)
105: PetscBLASIntCast(ds->n*ds->n+2*ds->n,&lwork);
106: PetscBLASIntCast(2*ds->n*ds->n+5*ds->n+1+n1,&lrwork);
107: #else
108: PetscBLASIntCast(2*ds->n*ds->n+6*ds->n+1,&lwork);
109: #endif
110: DSAllocateWork_Private(ds,lwork,lrwork,liwork);
111: work = ds->work;
112: iwork = ds->iwork;
113: off = ds->l+ds->l*ld;
114: MatDenseGetArray(ds->omat[DS_MAT_A],&A);
115: MatDenseGetArray(ds->omat[DS_MAT_B],&B);
116: MatDenseGetArray(ds->omat[DS_MAT_Q],&Q);
117: #if defined(PETSC_USE_COMPLEX)
118: rr = ds->rwork;
119: rwork = ds->rwork + n1;
120: lrwork = ds->lrwork - n1;
121: PetscCallBLAS("LAPACKsygvd",LAPACKsygvd_(&itype,"V","U",&n1,A+off,&ld,B+off,&ld,rr,work,&lwork,rwork,&lrwork,iwork,&liwork,&info));
122: for (i=0;i<n1;i++) wr[ds->l+i] = rr[i];
123: #else
124: PetscCallBLAS("LAPACKsygvd",LAPACKsygvd_(&itype,"V","U",&n1,A+off,&ld,B+off,&ld,wr+ds->l,work,&lwork,iwork,&liwork,&info));
125: #endif
126: SlepcCheckLapackInfo("sygvd",info);
127: PetscArrayzero(Q+ds->l*ld,n1*ld);
128: for (i=ds->l;i<ds->n;i++) PetscArraycpy(Q+ds->l+i*ld,A+ds->l+i*ld,n1);
129: PetscArrayzero(B+ds->l*ld,n1*ld);
130: PetscArrayzero(A+ds->l*ld,n1*ld);
131: for (i=ds->l;i<ds->n;i++) {
132: if (wi) wi[i] = 0.0;
133: B[i+i*ld] = 1.0;
134: A[i+i*ld] = wr[i];
135: }
136: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
137: MatDenseRestoreArray(ds->omat[DS_MAT_B],&B);
138: MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q);
139: return 0;
140: }
142: #if !defined(PETSC_HAVE_MPIUNI)
143: PetscErrorCode DSSynchronize_GHEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
144: {
145: PetscScalar *A,*B,*Q;
146: PetscInt ld=ds->ld,l=ds->l,k;
147: PetscMPIInt n,rank,off=0,size,ldn;
149: k = 2*(ds->n-l)*ld;
150: if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
151: if (eigr) k += (ds->n-l);
152: DSAllocateWork_Private(ds,k,0,0);
153: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
154: PetscMPIIntCast(ds->n-l,&n);
155: PetscMPIIntCast(ld*(ds->n-l),&ldn);
156: MatDenseGetArray(ds->omat[DS_MAT_A],&A);
157: MatDenseGetArray(ds->omat[DS_MAT_B],&B);
158: if (ds->state>DS_STATE_RAW) MatDenseGetArray(ds->omat[DS_MAT_Q],&Q);
159: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
160: if (!rank) {
161: MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
162: MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
163: if (ds->state>DS_STATE_RAW) MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
164: if (eigr) MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
165: }
166: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
167: if (rank) {
168: MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
169: MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
170: if (ds->state>DS_STATE_RAW) MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
171: if (eigr) MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
172: }
173: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
174: MatDenseRestoreArray(ds->omat[DS_MAT_B],&B);
175: if (ds->state>DS_STATE_RAW) MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q);
176: return 0;
177: }
178: #endif
180: PetscErrorCode DSHermitian_GHEP(DS ds,DSMatType m,PetscBool *flg)
181: {
182: if (m==DS_MAT_A || m==DS_MAT_B) *flg = PETSC_TRUE;
183: else *flg = PETSC_FALSE;
184: return 0;
185: }
187: /*MC
188: DSGHEP - Dense Generalized Hermitian Eigenvalue Problem.
190: Level: beginner
192: Notes:
193: The problem is expressed as A*X = B*X*Lambda, where both A and B are
194: real symmetric (or complex Hermitian) and B is positive-definite. Lambda
195: is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
196: After solve, A is overwritten with Lambda, and B is overwritten with I.
198: No intermediate state is implemented, nor compact storage.
200: Used DS matrices:
201: + DS_MAT_A - first problem matrix
202: . DS_MAT_B - second problem matrix
203: - DS_MAT_Q - matrix of B-orthogonal eigenvectors, which is equal to X
205: Implemented methods:
206: . 0 - Divide and Conquer (_sygvd)
208: .seealso: DSCreate(), DSSetType(), DSType
209: M*/
210: SLEPC_EXTERN PetscErrorCode DSCreate_GHEP(DS ds)
211: {
212: ds->ops->allocate = DSAllocate_GHEP;
213: ds->ops->view = DSView_GHEP;
214: ds->ops->vectors = DSVectors_GHEP;
215: ds->ops->solve[0] = DSSolve_GHEP;
216: ds->ops->sort = DSSort_GHEP;
217: #if !defined(PETSC_HAVE_MPIUNI)
218: ds->ops->synchronize = DSSynchronize_GHEP;
219: #endif
220: ds->ops->hermitian = DSHermitian_GHEP;
221: return 0;
222: }