Actual source code: qeplin.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }