Actual source code: qeplin.c

slepc-3.18.3 2023-03-24
Report Typos and Errors
  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:   MatGetSize(ctx->M,&M,&N);
 45:   MatGetLocalSize(ctx->M,&m,&n);
 46:   MatCreateConstantDiagonal(PetscObjectComm((PetscObject)ctx->M),m,n,M,N,1.0,&Id);
 47:   if (a!=0.0 && b!=0.0) {
 48:     MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
 49:     MatScale(T,-a*ctx->dsfactor*ctx->sfactor);
 50:     MatShift(T,b);
 51:   } else {
 52:     if (a==0.0) { T = Id; scalt = b; }
 53:     else { T = ctx->C; scalt = -a*ctx->dsfactor*ctx->sfactor; }
 54:   }
 55:   MatCreateTile(-b*ctx->dsfactor,ctx->K,a,Id,-ctx->dsfactor*a,ctx->K,scalt,T,A);
 56:   MatDestroy(&Id);
 57:   if (a!=0.0 && b!=0.0) MatDestroy(&T);
 58:   return 0;
 59: }

 61: PetscErrorCode MatCreateExplicit_Linear_NB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
 62: {
 63:   PetscInt       M,N,m,n;
 64:   Mat            Id,T=NULL;
 65:   PetscReal      a=ctx->alpha,b=ctx->beta;
 66:   PetscScalar    scalt=1.0;

 68:   MatGetSize(ctx->M,&M,&N);
 69:   MatGetLocalSize(ctx->M,&m,&n);
 70:   MatCreateConstantDiagonal(PetscObjectComm((PetscObject)ctx->M),m,n,M,N,1.0,&Id);
 71:   if (a!=0.0 && b!=0.0) {
 72:     MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
 73:     MatScale(T,b*ctx->dsfactor*ctx->sfactor);
 74:     MatShift(T,a);
 75:   } else {
 76:     if (b==0.0) { T = Id; scalt = a; }
 77:     else { T = ctx->C; scalt = b*ctx->dsfactor*ctx->sfactor; }
 78:   }
 79:   MatCreateTile(scalt,T,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,b,Id,a*ctx->sfactor*ctx->sfactor*ctx->dsfactor,ctx->M,B);
 80:   MatDestroy(&Id);
 81:   if (a!=0.0 && b!=0.0) MatDestroy(&T);
 82:   return 0;
 83: }

 85: /* - - - S - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 87: PetscErrorCode MatCreateExplicit_Linear_SA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
 88: {
 89:   Mat            T=NULL;
 90:   PetscScalar    scalt=1.0;
 91:   PetscReal      a=ctx->alpha,b=ctx->beta;

 93:   if (a!=0.0 && b!=0.0) {
 94:     MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
 95:     MatScale(T,a*ctx->dsfactor*ctx->sfactor);
 96:     MatAXPY(T,-b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,UNKNOWN_NONZERO_PATTERN);
 97:   } else {
 98:     if (a==0.0) { T = ctx->M; scalt = -b*ctx->dsfactor*ctx->sfactor*ctx->sfactor; }
 99:     else { T = ctx->C; scalt = a*ctx->dsfactor*ctx->sfactor; }
100:   }
101:   MatCreateTile(b*ctx->dsfactor,ctx->K,a*ctx->dsfactor,ctx->K,a*ctx->dsfactor,ctx->K,scalt,T,A);
102:   if (a!=0.0 && b!=0.0) MatDestroy(&T);
103:   return 0;
104: }

106: PetscErrorCode MatCreateExplicit_Linear_SB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
107: {
108:   Mat            T=NULL;
109:   PetscScalar    scalt=1.0;
110:   PetscReal      a=ctx->alpha,b=ctx->beta;

112:   if (a!=0.0 && b!=0.0) {
113:     MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
114:     MatScale(T,-b*ctx->dsfactor*ctx->sfactor);
115:     MatAXPY(T,a*ctx->dsfactor,ctx->K,UNKNOWN_NONZERO_PATTERN);
116:   } else {
117:     if (b==0.0) { T = ctx->K; scalt = a*ctx->dsfactor; }
118:     else { T = ctx->C; scalt = -b*ctx->dsfactor*ctx->sfactor; }
119:   }
120:   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);
121:   if (a!=0.0 && b!=0.0) MatDestroy(&T);
122:   return 0;
123: }

125: /* - - - H - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

127: PetscErrorCode MatCreateExplicit_Linear_HA(MPI_Comm comm,PEP_LINEAR *ctx,Mat *A)
128: {
129:   Mat            T=NULL;
130:   PetscScalar    scalt=1.0;
131:   PetscReal      a=ctx->alpha,b=ctx->beta;

133:   if (a!=0.0 && b!=0.0) {
134:     MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
135:     MatScale(T,a*ctx->dsfactor*ctx->sfactor);
136:     MatAXPY(T,b*ctx->dsfactor*ctx->sfactor*ctx->sfactor,ctx->M,UNKNOWN_NONZERO_PATTERN);
137:   } else {
138:     if (a==0.0) { T = ctx->M; scalt = b*ctx->dsfactor*ctx->sfactor*ctx->sfactor; }
139:     else { T = ctx->C; scalt = a*ctx->dsfactor*ctx->sfactor; }
140:   }
141:   MatCreateTile(a*ctx->dsfactor,ctx->K,-b*ctx->dsfactor,ctx->K,scalt,T,a*ctx->dsfactor,ctx->K,A);
142:   if (a!=0.0 && b!=0.0) MatDestroy(&T);
143:   return 0;
144: }

146: PetscErrorCode MatCreateExplicit_Linear_HB(MPI_Comm comm,PEP_LINEAR *ctx,Mat *B)
147: {
148:   Mat            T=NULL;
149:   PetscScalar    scalt=1.0;
150:   PetscReal      a=ctx->alpha,b=ctx->beta;

152:   if (a!=0.0 && b!=0.0) {
153:     MatDuplicate(ctx->C,MAT_COPY_VALUES,&T);
154:     MatScale(T,b*ctx->dsfactor*ctx->sfactor);
155:     MatAXPY(T,a*ctx->dsfactor,ctx->K,UNKNOWN_NONZERO_PATTERN);
156:   } else {
157:     if (b==0.0) { T = ctx->K; scalt = a*ctx->dsfactor; }
158:     else { T = ctx->C; scalt = b*ctx->dsfactor*ctx->sfactor; }
159:   }
160:   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);
161:   if (a!=0.0 && b!=0.0) MatDestroy(&T);
162:   return 0;
163: }