Actual source code: qeplin.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: Various types of linearization for quadratic eigenvalue problem
12: */
14: #include <slepc/private/pepimpl.h>
15: #include "linear.h"
17: /*
18: Given the quadratic problem (l^2*M + l*C + K)*x = 0 the linearization is
19: A*z = l*B*z for z = [ x ] and A,B defined as follows:
20: [ l*x ]
22: N:
23: A = [-bK aI ] B = [ aI+bC bM ]
24: [-aK -aC+bI ] [ bI aM ]
26: S:
27: A = [ bK aK ] B = [ aK-bC -bM ]
28: [ aK aC-bM ] [-bM -aM ]
30: H:
31: A = [ aK -bK ] B = [ bM aK+bC ]
32: [ aC+bM aK ] [-aM bM ]
33: */
35: /* - - - N - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37: PetscErrorCode MatCreateExplicit_Linear_NA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
38: {
39: PetscInt M,N,m,n;
40: Mat Id,T=NULL;
41: PetscReal a=ctx->alpha,b=ctx->beta;
42: PetscScalar scalt=1.0;
44: PetscFunctionBegin;
45: PetscCall(MatGetSize(ctx->M,&M,&N));
46: PetscCall(MatGetLocalSize(ctx->M,&m,&n));
47: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)ctx->M),m,n,M,N,1.0,&Id));
48: if (a!=0.0 && b!=0.0) {
49: PetscCall(MatDuplicate(ctx->C,MAT_COPY_VALUES,&T));
50: PetscCall(MatScale(T,-a*ctx->dsfactor*ctx->sfactor));
51: PetscCall(MatShift(T,b));
52: } else {
53: if (a==0.0) { T = Id; scalt = b; }
54: else { T = ctx->C; scalt = -a*ctx->dsfactor*ctx->sfactor; }
55: }
56: PetscCall(MatCreateTile(-b*ctx->dsfactor,ctx->K,a,Id,-ctx->dsfactor*a,ctx->K,scalt,T,A));
57: PetscCall(MatDestroy(&Id));
58: if (a!=0.0 && b!=0.0) PetscCall(MatDestroy(&T));
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: PetscErrorCode MatCreateExplicit_Linear_NB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
63: {
64: PetscInt M,N,m,n;
65: Mat Id,T=NULL;
66: PetscReal a=ctx->alpha,b=ctx->beta;
67: PetscScalar scalt=1.0;
69: PetscFunctionBegin;
70: PetscCall(MatGetSize(ctx->M,&M,&N));
71: PetscCall(MatGetLocalSize(ctx->M,&m,&n));
72: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)ctx->M),m,n,M,N,1.0,&Id));
73: if (a!=0.0 && b!=0.0) {
74: PetscCall(MatDuplicate(ctx->C,MAT_COPY_VALUES,&T));
75: PetscCall(MatScale(T,b*ctx->dsfactor*ctx->sfactor));
76: PetscCall(MatShift(T,a));
77: } else {
78: if (b==0.0) { T = Id; scalt = a; }
79: else { T = ctx->C; scalt = b*ctx->dsfactor*ctx->sfactor; }
80: }
81: PetscCall(MatCreateTile(scalt,T,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,b,Id,a*ctx->sfactor*ctx->sfactor*ctx->dsfactor,ctx->M,B));
82: PetscCall(MatDestroy(&Id));
83: if (a!=0.0 && b!=0.0) PetscCall(MatDestroy(&T));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /* - - - S - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: PetscErrorCode MatCreateExplicit_Linear_SA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
90: {
91: Mat T=NULL;
92: PetscScalar scalt=1.0;
93: PetscReal a=ctx->alpha,b=ctx->beta;
95: PetscFunctionBegin;
96: if (a!=0.0 && b!=0.0) {
97: PetscCall(MatDuplicate(ctx->C,MAT_COPY_VALUES,&T));
98: PetscCall(MatScale(T,a*ctx->dsfactor*ctx->sfactor));
99: PetscCall(MatAXPY(T,-b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,UNKNOWN_NONZERO_PATTERN));
100: } else {
101: if (a==0.0) { T = ctx->M; scalt = -b*ctx->dsfactor*ctx->sfactor*ctx->sfactor; }
102: else { T = ctx->C; scalt = a*ctx->dsfactor*ctx->sfactor; }
103: }
104: PetscCall(MatCreateTile(b*ctx->dsfactor,ctx->K,a*ctx->dsfactor,ctx->K,a*ctx->dsfactor,ctx->K,scalt,T,A));
105: if (a!=0.0 && b!=0.0) PetscCall(MatDestroy(&T));
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: PetscErrorCode MatCreateExplicit_Linear_SB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
110: {
111: Mat T=NULL;
112: PetscScalar scalt=1.0;
113: PetscReal a=ctx->alpha,b=ctx->beta;
115: PetscFunctionBegin;
116: if (a!=0.0 && b!=0.0) {
117: PetscCall(MatDuplicate(ctx->C,MAT_COPY_VALUES,&T));
118: PetscCall(MatScale(T,-b*ctx->dsfactor*ctx->sfactor));
119: PetscCall(MatAXPY(T,a*ctx->dsfactor,ctx->K,UNKNOWN_NONZERO_PATTERN));
120: } else {
121: if (b==0.0) { T = ctx->K; scalt = a*ctx->dsfactor; }
122: else { T = ctx->C; scalt = -b*ctx->dsfactor*ctx->sfactor; }
123: }
124: PetscCall(MatCreateTile(scalt,T,-b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,-b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,-a*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,B));
125: if (a!=0.0 && b!=0.0) PetscCall(MatDestroy(&T));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /* - - - H - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: PetscErrorCode MatCreateExplicit_Linear_HA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
132: {
133: Mat T=NULL;
134: PetscScalar scalt=1.0;
135: PetscReal a=ctx->alpha,b=ctx->beta;
137: PetscFunctionBegin;
138: if (a!=0.0 && b!=0.0) {
139: PetscCall(MatDuplicate(ctx->C,MAT_COPY_VALUES,&T));
140: PetscCall(MatScale(T,a*ctx->dsfactor*ctx->sfactor));
141: PetscCall(MatAXPY(T,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,UNKNOWN_NONZERO_PATTERN));
142: } else {
143: if (a==0.0) { T = ctx->M; scalt = b*ctx->dsfactor*ctx->sfactor*ctx->sfactor; }
144: else { T = ctx->C; scalt = a*ctx->dsfactor*ctx->sfactor; }
145: }
146: PetscCall(MatCreateTile(a*ctx->dsfactor,ctx->K,-b*ctx->dsfactor,ctx->K,scalt,T,a*ctx->dsfactor,ctx->K,A));
147: if (a!=0.0 && b!=0.0) PetscCall(MatDestroy(&T));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: PetscErrorCode MatCreateExplicit_Linear_HB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
152: {
153: Mat T=NULL;
154: PetscScalar scalt=1.0;
155: PetscReal a=ctx->alpha,b=ctx->beta;
157: PetscFunctionBegin;
158: if (a!=0.0 && b!=0.0) {
159: PetscCall(MatDuplicate(ctx->C,MAT_COPY_VALUES,&T));
160: PetscCall(MatScale(T,b*ctx->dsfactor*ctx->sfactor));
161: PetscCall(MatAXPY(T,a*ctx->dsfactor,ctx->K,UNKNOWN_NONZERO_PATTERN));
162: } else {
163: if (b==0.0) { T = ctx->K; scalt = a*ctx->dsfactor; }
164: else { T = ctx->C; scalt = b*ctx->dsfactor*ctx->sfactor; }
165: }
166: PetscCall(MatCreateTile(b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,scalt,T,-a*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,B));
167: if (a!=0.0 && b!=0.0) PetscCall(MatDestroy(&T));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }