Actual source code: qeplin.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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 "linearp.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 ]
27: S:
28: A = [ bK aK ] B = [ aK-bC -bM ]
29: [ aK aC-bM ] [-bM -aM ]
32: H:
33: A = [ aK -bK ] B = [ bM aK+bC ]
34: [ aC+bM aK ] [-aM bM ]
35: */
37: /* - - - N - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39: PetscErrorCode MatCreateExplicit_Linear_NA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
40: {
42: PetscInt M,N,m,n,i,Istart,Iend;
43: Mat Id,T=NULL;
44: PetscReal a=ctx->alpha,b=ctx->beta;
45: PetscScalar scalt=1.0;
48: MatGetSize(ctx->M,&M,&N);
49: MatGetLocalSize(ctx->M,&m,&n);
50: MatCreate(PetscObjectComm((PetscObject)ctx->M),&Id);
51: MatSetSizes(Id,m,n,M,N);
52: MatSetFromOptions(Id);
53: MatSetUp(Id);
54: MatGetOwnershipRange(Id,&Istart,&Iend);
55: for (i=Istart;i<Iend;i++) {
56: MatSetValue(Id,i,i,1.0,INSERT_VALUES);
57: }
58: MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY);
60: if (a!=0.0 && b!=0.0) {
61: MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
62: MatScale(T,-a*ctx->dsfactor*ctx->sfactor);
63: MatShift(T,b);
64: } else {
65: if (a==0.0) { T = Id; scalt = b; }
66: else { T = ctx->C; scalt = -a*ctx->dsfactor*ctx->sfactor; }
67: }
68: MatCreateTile(-b*ctx->dsfactor,ctx->K,a,Id,-ctx->dsfactor*a,ctx->K,scalt,T,A);
69: MatDestroy(&Id);
70: if (a!=0.0 && b!=0.0) {
71: MatDestroy(&T);
72: }
73: return(0);
74: }
76: PetscErrorCode MatCreateExplicit_Linear_NB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
77: {
79: PetscInt M,N,m,n,i,Istart,Iend;
80: Mat Id,T=NULL;
81: PetscReal a=ctx->alpha,b=ctx->beta;
82: PetscScalar scalt=1.0;
85: MatGetSize(ctx->M,&M,&N);
86: MatGetLocalSize(ctx->M,&m,&n);
87: MatCreate(PetscObjectComm((PetscObject)ctx->M),&Id);
88: MatSetSizes(Id,m,n,M,N);
89: MatSetFromOptions(Id);
90: MatSetUp(Id);
91: MatGetOwnershipRange(Id,&Istart,&Iend);
92: for (i=Istart;i<Iend;i++) {
93: MatSetValue(Id,i,i,1.0,INSERT_VALUES);
94: }
95: MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY);
97: if (a!=0.0 && b!=0.0) {
98: MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
99: MatScale(T,b*ctx->dsfactor*ctx->sfactor);
100: MatShift(T,a);
101: } else {
102: if (b==0.0) { T = Id; scalt = a; }
103: else { T = ctx->C; scalt = b*ctx->dsfactor*ctx->sfactor; }
104: }
105: MatCreateTile(scalt,T,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,b,Id,a*ctx->sfactor*ctx->sfactor*ctx->dsfactor,ctx->M,B);
106: MatDestroy(&Id);
107: if (a!=0.0 && b!=0.0) {
108: MatDestroy(&T);
109: }
110: return(0);
111: }
113: /* - - - S - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: PetscErrorCode MatCreateExplicit_Linear_SA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
116: {
118: Mat T=NULL;
119: PetscScalar scalt=1.0;
120: PetscReal a=ctx->alpha,b=ctx->beta;
123: if (a!=0.0 && b!=0.0) {
124: MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
125: MatScale(T,a*ctx->dsfactor*ctx->sfactor);
126: MatAXPY(T,-b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,DIFFERENT_NONZERO_PATTERN);
127: } else {
128: if (a==0.0) { T = ctx->M; scalt = -b*ctx->dsfactor*ctx->sfactor*ctx->sfactor; }
129: else { T = ctx->C; scalt = a*ctx->dsfactor*ctx->sfactor; }
130: }
131: MatCreateTile(b*ctx->dsfactor,ctx->K,a*ctx->dsfactor,ctx->K,a*ctx->dsfactor,ctx->K,scalt,T,A);
132: if (a!=0.0 && b!=0.0) {
133: MatDestroy(&T);
134: }
135: return(0);
136: }
138: PetscErrorCode MatCreateExplicit_Linear_SB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
139: {
141: Mat T=NULL;
142: PetscScalar scalt=1.0;
143: PetscReal a=ctx->alpha,b=ctx->beta;
146: if (a!=0.0 && b!=0.0) {
147: MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
148: MatScale(T,-b*ctx->dsfactor*ctx->sfactor);
149: MatAXPY(T,a*ctx->dsfactor,ctx->K,DIFFERENT_NONZERO_PATTERN);
150: } else {
151: if (b==0.0) { T = ctx->K; scalt = a*ctx->dsfactor; }
152: else { T = ctx->C; scalt = -b*ctx->dsfactor*ctx->sfactor; }
153: }
154: 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);
155: if (a!=0.0 && b!=0.0) {
156: MatDestroy(&T);
157: }
158: return(0);
159: }
161: /* - - - H - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: PetscErrorCode MatCreateExplicit_Linear_HA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
164: {
166: Mat T=NULL;
167: PetscScalar scalt=1.0;
168: PetscReal a=ctx->alpha,b=ctx->beta;
171: if (a!=0.0 && b!=0.0) {
172: MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
173: MatScale(T,a*ctx->dsfactor*ctx->sfactor);
174: MatAXPY(T,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,DIFFERENT_NONZERO_PATTERN);
175: } else {
176: if (a==0.0) { T = ctx->M; scalt = b*ctx->dsfactor*ctx->sfactor*ctx->sfactor; }
177: else { T = ctx->C; scalt = a*ctx->dsfactor*ctx->sfactor; }
178: }
179: MatCreateTile(a*ctx->dsfactor,ctx->K,-b*ctx->dsfactor,ctx->K,scalt,T,a*ctx->dsfactor,ctx->K,A);
180: if (a!=0.0 && b!=0.0) {
181: MatDestroy(&T);
182: }
183: return(0);
184: }
186: PetscErrorCode MatCreateExplicit_Linear_HB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
187: {
189: Mat T=NULL;
190: PetscScalar scalt=1.0;
191: PetscReal a=ctx->alpha,b=ctx->beta;
194: if (a!=0.0 && b!=0.0) {
195: MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
196: MatScale(T,b*ctx->dsfactor*ctx->sfactor);
197: MatAXPY(T,a*ctx->dsfactor,ctx->K,DIFFERENT_NONZERO_PATTERN);
198: } else {
199: if (b==0.0) { T = ctx->K; scalt = a*ctx->dsfactor; }
200: else { T = ctx->C; scalt = b*ctx->dsfactor*ctx->sfactor; }
201: }
202: 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);
203: if (a!=0.0 && b!=0.0) {
204: MatDestroy(&T);
205: }
206: return(0);
207: }