Actual source code: peputils.c
slepc-3.21.1 2024-04-26
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: */
10: /*
11: Common subroutines for several PEP solvers
12: */
14: #include <slepc/private/pepimpl.h>
15: #include <slepcblaslapack.h>
17: /*
18: Computes T_j=phy_idx(T) for a matrix T.
19: Tp (if j>0) and Tpp (if j>1) are the evaluation
20: of phy_(j-1) and phy(j-2)respectively.
21: */
22: PetscErrorCode PEPEvaluateBasisMat(PEP pep,PetscInt k,PetscScalar *T,PetscInt ldt,PetscInt idx,PetscScalar *Tpp,PetscInt ldtpp,PetscScalar *Tp,PetscInt ldtp,PetscScalar *Tj,PetscInt ldtj)
23: {
24: PetscInt i;
25: PetscReal *ca,*cb,*cg;
26: PetscScalar g,a;
27: PetscBLASInt k_,ldt_,ldtj_,ldtp_;
29: PetscFunctionBegin;
30: if (idx==0) {
31: for (i=0;i<k;i++) {
32: PetscCall(PetscArrayzero(Tj+i*ldtj,k));
33: Tj[i+i*ldtj] = 1.0;
34: }
35: } else {
36: PetscCall(PetscBLASIntCast(ldt,&ldt_));
37: PetscCall(PetscBLASIntCast(ldtj,&ldtj_));
38: PetscCall(PetscBLASIntCast(ldtp,&ldtp_));
39: PetscCall(PetscBLASIntCast(k,&k_));
40: ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
41: for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
42: if (idx>1) {
43: for (i=0;i<k;i++) PetscCall(PetscArraycpy(Tj+i*ldtj,Tpp+i*ldtpp,k));
44: }
45: a = 1/ca[idx-1];
46: g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
47: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,Tp,&ldtp_,&g,Tj,&ldtj_));
48: for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*
54: PEPEvaluateBasis - evaluate the polynomial basis on a given parameter sigma
55: */
56: PetscErrorCode PEPEvaluateBasis(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
57: {
58: PetscInt nmat=pep->nmat,k;
59: PetscReal *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
61: PetscFunctionBegin;
62: if (ivals) for (k=0;k<nmat;k++) ivals[k] = 0.0;
63: vals[0] = 1.0;
64: vals[1] = (sigma-b[0])/a[0];
65: #if !defined(PETSC_USE_COMPLEX)
66: if (ivals) ivals[1] = isigma/a[0];
67: #endif
68: for (k=2;k<nmat;k++) {
69: vals[k] = ((sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2])/a[k-1];
70: if (ivals) vals[k] -= isigma*ivals[k-1]/a[k-1];
71: #if !defined(PETSC_USE_COMPLEX)
72: if (ivals) ivals[k] = ((sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2])/a[k-1];
73: #endif
74: }
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*
79: PEPEvaluateBasisDerivative - evaluate the derivative of the polynomial basis on a given parameter sigma
80: */
81: PetscErrorCode PEPEvaluateBasisDerivative(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
82: {
83: PetscInt nmat=pep->nmat,k;
84: PetscReal *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
86: PetscFunctionBegin;
87: PetscCall(PEPEvaluateBasis(pep,sigma,isigma,vals,ivals));
88: for (k=nmat-1;k>0;k--) {
89: vals[k] = vals[k-1];
90: if (ivals) ivals[k] = ivals[k-1];
91: }
92: vals[0] = 0.0;
93: vals[1] = vals[1]/a[0];
94: #if !defined(PETSC_USE_COMPLEX)
95: if (ivals) ivals[1] = ivals[1]/a[0];
96: #endif
97: for (k=2;k<nmat;k++) {
98: vals[k] += (sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2];
99: if (ivals) vals[k] -= isigma*ivals[k-1];
100: vals[k] /= a[k-1];
101: #if !defined(PETSC_USE_COMPLEX)
102: if (ivals) {
103: ivals[k] += (sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2];
104: ivals[k] /= a[k-1];
105: }
106: #endif
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }