| 1910 |
jroman |
1 |
/*
|
|
|
2 |
|
|
|
3 |
Linearization for general QEP, companion form 2.
|
|
|
4 |
|
|
|
5 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
6 |
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
| 2116 |
eromero |
7 |
Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
|
| 1910 |
jroman |
8 |
|
|
|
9 |
This file is part of SLEPc.
|
|
|
10 |
|
|
|
11 |
SLEPc is free software: you can redistribute it and/or modify it under the
|
|
|
12 |
terms of version 3 of the GNU Lesser General Public License as published by
|
|
|
13 |
the Free Software Foundation.
|
|
|
14 |
|
|
|
15 |
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
|
16 |
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
|
17 |
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
|
|
18 |
more details.
|
|
|
19 |
|
|
|
20 |
You should have received a copy of the GNU Lesser General Public License
|
|
|
21 |
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
|
|
|
22 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
23 |
*/
|
|
|
24 |
|
|
|
25 |
#include "private/qepimpl.h" /*I "slepcqep.h" I*/
|
|
|
26 |
#include "slepceps.h"
|
|
|
27 |
#include "linearp.h"
|
|
|
28 |
|
|
|
29 |
/*
|
|
|
30 |
Given the quadratic problem (l^2*M + l*C + K)*x = 0 the following
|
|
|
31 |
linearization is employed:
|
|
|
32 |
|
|
|
33 |
A*z = l*B*z where A = [ -K 0 ] B = [ C M ] z = [ x ]
|
|
|
34 |
[ 0 I ] [ I 0 ] [ l*x ]
|
|
|
35 |
*/
|
|
|
36 |
|
|
|
37 |
#undef __FUNCT__
|
|
|
38 |
#define __FUNCT__ "MatMult_QEPLINEAR_N2A"
|
|
|
39 |
PetscErrorCode MatMult_QEPLINEAR_N2A(Mat A,Vec x,Vec y)
|
|
|
40 |
{
|
|
|
41 |
PetscErrorCode ierr;
|
|
|
42 |
QEP_LINEAR *ctx;
|
|
|
43 |
PetscScalar *px,*py;
|
|
|
44 |
PetscInt m;
|
|
|
45 |
|
|
|
46 |
PetscFunctionBegin;
|
|
|
47 |
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
|
|
|
48 |
ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr);
|
|
|
49 |
ierr = VecGetArray(x,&px);CHKERRQ(ierr);
|
|
|
50 |
ierr = VecGetArray(y,&py);CHKERRQ(ierr);
|
|
|
51 |
ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr);
|
|
|
52 |
ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr);
|
|
|
53 |
ierr = VecPlaceArray(ctx->y1,py);CHKERRQ(ierr);
|
|
|
54 |
ierr = VecPlaceArray(ctx->y2,py+m);CHKERRQ(ierr);
|
|
|
55 |
/* y1 = -K*x1 */
|
|
|
56 |
ierr = MatMult(ctx->K,ctx->x1,ctx->y1);CHKERRQ(ierr);
|
|
|
57 |
ierr = VecScale(ctx->y1,-1.0);CHKERRQ(ierr);
|
|
|
58 |
/* y2 = x2 */
|
|
|
59 |
ierr = VecCopy(ctx->x2,ctx->y2);CHKERRQ(ierr);
|
|
|
60 |
ierr = VecResetArray(ctx->x1);CHKERRQ(ierr);
|
|
|
61 |
ierr = VecResetArray(ctx->x2);CHKERRQ(ierr);
|
|
|
62 |
ierr = VecResetArray(ctx->y1);CHKERRQ(ierr);
|
|
|
63 |
ierr = VecResetArray(ctx->y2);CHKERRQ(ierr);
|
|
|
64 |
ierr = VecRestoreArray(x,&px);CHKERRQ(ierr);
|
|
|
65 |
ierr = VecRestoreArray(y,&py);CHKERRQ(ierr);
|
|
|
66 |
PetscFunctionReturn(0);
|
|
|
67 |
}
|
|
|
68 |
|
|
|
69 |
#undef __FUNCT__
|
|
|
70 |
#define __FUNCT__ "MatMult_QEPLINEAR_N2B"
|
|
|
71 |
PetscErrorCode MatMult_QEPLINEAR_N2B(Mat B,Vec x,Vec y)
|
|
|
72 |
{
|
|
|
73 |
PetscErrorCode ierr;
|
|
|
74 |
QEP_LINEAR *ctx;
|
|
|
75 |
PetscScalar *px,*py;
|
|
|
76 |
PetscInt m;
|
|
|
77 |
|
|
|
78 |
PetscFunctionBegin;
|
|
|
79 |
ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr);
|
|
|
80 |
ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr);
|
|
|
81 |
ierr = VecGetArray(x,&px);CHKERRQ(ierr);
|
|
|
82 |
ierr = VecGetArray(y,&py);CHKERRQ(ierr);
|
|
|
83 |
ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr);
|
|
|
84 |
ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr);
|
|
|
85 |
ierr = VecPlaceArray(ctx->y1,py);CHKERRQ(ierr);
|
|
|
86 |
ierr = VecPlaceArray(ctx->y2,py+m);CHKERRQ(ierr);
|
|
|
87 |
/* y1 = C*x1 + M*x2 */
|
|
|
88 |
ierr = MatMult(ctx->C,ctx->x1,ctx->y1);CHKERRQ(ierr);
|
| 1919 |
jroman |
89 |
ierr = VecScale(ctx->y1,ctx->sfactor);CHKERRQ(ierr);
|
| 1910 |
jroman |
90 |
ierr = MatMult(ctx->M,ctx->x2,ctx->y2);CHKERRQ(ierr);
|
| 1919 |
jroman |
91 |
ierr = VecAXPY(ctx->y1,ctx->sfactor*ctx->sfactor,ctx->y2);CHKERRQ(ierr);
|
| 1910 |
jroman |
92 |
/* y2 = x1 */
|
|
|
93 |
ierr = VecCopy(ctx->x1,ctx->y2);CHKERRQ(ierr);
|
|
|
94 |
ierr = VecResetArray(ctx->x1);CHKERRQ(ierr);
|
|
|
95 |
ierr = VecResetArray(ctx->x2);CHKERRQ(ierr);
|
|
|
96 |
ierr = VecResetArray(ctx->y1);CHKERRQ(ierr);
|
|
|
97 |
ierr = VecResetArray(ctx->y2);CHKERRQ(ierr);
|
|
|
98 |
ierr = VecRestoreArray(x,&px);CHKERRQ(ierr);
|
|
|
99 |
ierr = VecRestoreArray(y,&py);CHKERRQ(ierr);
|
|
|
100 |
PetscFunctionReturn(0);
|
|
|
101 |
}
|
|
|
102 |
|
|
|
103 |
#undef __FUNCT__
|
|
|
104 |
#define __FUNCT__ "MatGetDiagonal_QEPLINEAR_N2A"
|
|
|
105 |
PetscErrorCode MatGetDiagonal_QEPLINEAR_N2A(Mat A,Vec diag)
|
|
|
106 |
{
|
|
|
107 |
PetscErrorCode ierr;
|
|
|
108 |
QEP_LINEAR *ctx;
|
|
|
109 |
PetscScalar *pd;
|
|
|
110 |
PetscInt m;
|
|
|
111 |
|
|
|
112 |
PetscFunctionBegin;
|
|
|
113 |
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
|
|
|
114 |
ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr);
|
|
|
115 |
ierr = VecGetArray(diag,&pd);CHKERRQ(ierr);
|
|
|
116 |
ierr = VecPlaceArray(ctx->x1,pd);CHKERRQ(ierr);
|
|
|
117 |
ierr = VecPlaceArray(ctx->x2,pd+m);CHKERRQ(ierr);
|
|
|
118 |
ierr = MatGetDiagonal(ctx->K,ctx->x1);CHKERRQ(ierr);
|
|
|
119 |
ierr = VecScale(ctx->x1,-1.0);CHKERRQ(ierr);
|
|
|
120 |
ierr = VecSet(ctx->x2,1.0);CHKERRQ(ierr);
|
|
|
121 |
ierr = VecResetArray(ctx->x1);CHKERRQ(ierr);
|
|
|
122 |
ierr = VecResetArray(ctx->x2);CHKERRQ(ierr);
|
|
|
123 |
ierr = VecRestoreArray(diag,&pd);CHKERRQ(ierr);
|
|
|
124 |
PetscFunctionReturn(0);
|
|
|
125 |
}
|
|
|
126 |
|
|
|
127 |
#undef __FUNCT__
|
|
|
128 |
#define __FUNCT__ "MatGetDiagonal_QEPLINEAR_N2B"
|
|
|
129 |
PetscErrorCode MatGetDiagonal_QEPLINEAR_N2B(Mat B,Vec diag)
|
|
|
130 |
{
|
|
|
131 |
PetscErrorCode ierr;
|
|
|
132 |
QEP_LINEAR *ctx;
|
|
|
133 |
PetscScalar *pd;
|
|
|
134 |
PetscInt m;
|
|
|
135 |
|
|
|
136 |
PetscFunctionBegin;
|
|
|
137 |
ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr);
|
|
|
138 |
ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr);
|
|
|
139 |
ierr = VecGetArray(diag,&pd);CHKERRQ(ierr);
|
|
|
140 |
ierr = VecPlaceArray(ctx->x1,pd);CHKERRQ(ierr);
|
|
|
141 |
ierr = VecPlaceArray(ctx->x2,pd+m);CHKERRQ(ierr);
|
|
|
142 |
ierr = MatGetDiagonal(ctx->C,ctx->x1);CHKERRQ(ierr);
|
| 1919 |
jroman |
143 |
ierr = VecScale(ctx->x1,ctx->sfactor);CHKERRQ(ierr);
|
| 1910 |
jroman |
144 |
ierr = VecSet(ctx->x2,0.0);CHKERRQ(ierr);
|
|
|
145 |
ierr = VecResetArray(ctx->x1);CHKERRQ(ierr);
|
|
|
146 |
ierr = VecResetArray(ctx->x2);CHKERRQ(ierr);
|
|
|
147 |
ierr = VecRestoreArray(diag,&pd);CHKERRQ(ierr);
|
|
|
148 |
PetscFunctionReturn(0);
|
|
|
149 |
}
|
|
|
150 |
|
|
|
151 |
#undef __FUNCT__
|
|
|
152 |
#define __FUNCT__ "MatCreateExplicit_QEPLINEAR_N2A"
|
|
|
153 |
PetscErrorCode MatCreateExplicit_QEPLINEAR_N2A(MPI_Comm comm,QEP_LINEAR *ctx,Mat *A)
|
|
|
154 |
{
|
|
|
155 |
PetscErrorCode ierr;
|
|
|
156 |
PetscInt M,N,m,n,i,start,end,ncols;
|
|
|
157 |
const PetscInt *cols;
|
|
|
158 |
const PetscScalar *vals;
|
|
|
159 |
|
|
|
160 |
PetscFunctionBegin;
|
|
|
161 |
ierr = MatGetSize(ctx->M,&M,&N);CHKERRQ(ierr);
|
|
|
162 |
ierr = MatGetLocalSize(ctx->M,&m,&n);CHKERRQ(ierr);
|
|
|
163 |
ierr = MatCreate(comm,A);CHKERRQ(ierr);
|
|
|
164 |
ierr = MatSetSizes(*A,m+n,m+n,M+N,M+N);CHKERRQ(ierr);
|
|
|
165 |
ierr = MatSetFromOptions(*A);CHKERRQ(ierr);
|
|
|
166 |
ierr = MatGetOwnershipRange(ctx->M,&start,&end);CHKERRQ(ierr);
|
|
|
167 |
for (i=start;i<end;i++) {
|
|
|
168 |
ierr = MatGetRow(ctx->K,i,&ncols,&cols,&vals);CHKERRQ(ierr);
|
|
|
169 |
ierr = MatSetValues(*A,1,&i,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
|
|
|
170 |
ierr = MatRestoreRow(ctx->K,i,&ncols,&cols,&vals);CHKERRQ(ierr);
|
|
|
171 |
ierr = MatSetValue(*A,i+M,i+M,-1.0,INSERT_VALUES);CHKERRQ(ierr);
|
|
|
172 |
}
|
|
|
173 |
ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
|
|
|
174 |
ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
|
|
|
175 |
ierr = MatScale(*A,-1.0);CHKERRQ(ierr);
|
|
|
176 |
PetscFunctionReturn(0);
|
|
|
177 |
}
|
|
|
178 |
|
|
|
179 |
#undef __FUNCT__
|
|
|
180 |
#define __FUNCT__ "MatCreateExplicit_QEPLINEAR_N2B"
|
|
|
181 |
PetscErrorCode MatCreateExplicit_QEPLINEAR_N2B(MPI_Comm comm,QEP_LINEAR *ctx,Mat *B)
|
|
|
182 |
{
|
|
|
183 |
PetscErrorCode ierr;
|
|
|
184 |
PetscInt M,N,m,n,i,j,start,end,ncols,*pos;
|
| 1919 |
jroman |
185 |
PetscScalar *svals;
|
| 1910 |
jroman |
186 |
const PetscInt *cols;
|
|
|
187 |
const PetscScalar *vals;
|
|
|
188 |
|
|
|
189 |
PetscFunctionBegin;
|
|
|
190 |
ierr = MatGetSize(ctx->M,&M,&N);CHKERRQ(ierr);
|
|
|
191 |
ierr = MatGetLocalSize(ctx->M,&m,&n);CHKERRQ(ierr);
|
|
|
192 |
ierr = MatCreate(comm,B);CHKERRQ(ierr);
|
|
|
193 |
ierr = MatSetSizes(*B,m+n,m+n,M+N,M+N);CHKERRQ(ierr);
|
|
|
194 |
ierr = MatSetFromOptions(*B);CHKERRQ(ierr);
|
|
|
195 |
ierr = PetscMalloc(sizeof(PetscInt)*n,&pos);CHKERRQ(ierr);
|
| 1919 |
jroman |
196 |
ierr = PetscMalloc(sizeof(PetscScalar)*n,&svals);CHKERRQ(ierr);
|
| 1910 |
jroman |
197 |
ierr = MatGetOwnershipRange(ctx->M,&start,&end);CHKERRQ(ierr);
|
|
|
198 |
for (i=start;i<end;i++) {
|
|
|
199 |
ierr = MatGetRow(ctx->C,i,&ncols,&cols,&vals);CHKERRQ(ierr);
|
| 1919 |
jroman |
200 |
for (j=0;j<ncols;j++) {
|
|
|
201 |
svals[j] = vals[j]*ctx->sfactor;
|
|
|
202 |
}
|
|
|
203 |
ierr = MatSetValues(*B,1,&i,ncols,cols,svals,INSERT_VALUES);CHKERRQ(ierr);
|
| 1910 |
jroman |
204 |
ierr = MatRestoreRow(ctx->C,i,&ncols,&cols,&vals);CHKERRQ(ierr);
|
|
|
205 |
ierr = MatGetRow(ctx->M,i,&ncols,&cols,&vals);CHKERRQ(ierr);
|
| 1919 |
jroman |
206 |
for (j=0;j<ncols;j++) {
|
| 1910 |
jroman |
207 |
pos[j] = cols[j] + M;
|
| 1919 |
jroman |
208 |
svals[j] = vals[j]*ctx->sfactor*ctx->sfactor;
|
|
|
209 |
}
|
|
|
210 |
ierr = MatSetValues(*B,1,&i,ncols,pos,svals,INSERT_VALUES);CHKERRQ(ierr);
|
| 1910 |
jroman |
211 |
ierr = MatRestoreRow(ctx->M,i,&ncols,&cols,&vals);CHKERRQ(ierr);
|
|
|
212 |
ierr = MatSetValue(*B,i+M,i,1.0,INSERT_VALUES);CHKERRQ(ierr);
|
|
|
213 |
}
|
|
|
214 |
ierr = PetscFree(pos);CHKERRQ(ierr);
|
| 1919 |
jroman |
215 |
ierr = PetscFree(svals);CHKERRQ(ierr);
|
| 1910 |
jroman |
216 |
ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
|
|
|
217 |
ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
|
|
|
218 |
PetscFunctionReturn(0);
|
|
|
219 |
}
|
|
|
220 |
|