WebSVN - slepc-dev - Blame - Rev 2116 - /trunk/src/st/impls/fold/fold.c

Rev	Author	Line No.	Line
778	dsic.upv.es!antodo	1	/*
804	dsic.upv.es!antodo	2	Folding spectral transformation, applies (A + sigma I)^2 as operator, or
		3	inv(B)(A + sigma I)^2 for generalized problems
1376	slepc	4
		5	- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1672	slepc	6	SLEPc - Scalable Library for Eigenvalue Problem Computations
2116	eromero	7	Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
1376	slepc	8
1672	slepc	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/>.
1376	slepc	22	- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
778	dsic.upv.es!antodo	23	*/
1376	slepc	24
1521	slepc	25	#include "private/stimpl.h" /I "slepcst.h" I/
778	dsic.upv.es!antodo	26
		27	typedef struct {
		28	Vec w2;
		29	} ST_FOLD;
		30
		31	#undef __FUNCT__
		32	#define __FUNCT__ "STApply_Fold"
		33	PetscErrorCode STApply_Fold(ST st,Vec x,Vec y)
		34	{
		35	PetscErrorCode ierr;
		36	ST_FOLD ctx = (ST_FOLD ) st->data;
		37
		38	PetscFunctionBegin;
		39	if (st->B) {
804	dsic.upv.es!antodo	40	/* generalized eigenproblem: y = (B^-1 A + sI)^2 x */
828	dsic.upv.es!antodo	41	ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
		42	ierr = STAssociatedKSPSolve(st,st->w,ctx->w2);CHKERRQ(ierr);
778	dsic.upv.es!antodo	43	if (st->sigma != 0.0) {
828	dsic.upv.es!antodo	44	ierr = VecAXPY(ctx->w2,-st->sigma,x);CHKERRQ(ierr);
778	dsic.upv.es!antodo	45	}
828	dsic.upv.es!antodo	46	ierr = MatMult(st->A,ctx->w2,st->w);CHKERRQ(ierr);
		47	ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
		48	if (st->sigma != 0.0) {
		49	ierr = VecAXPY(y,-st->sigma,ctx->w2);CHKERRQ(ierr);
		50	}
778	dsic.upv.es!antodo	51	} else {
804	dsic.upv.es!antodo	52	/* standard eigenproblem: y = (A + sI)^2 x */
778	dsic.upv.es!antodo	53	ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
		54	if (st->sigma != 0.0) {
828	dsic.upv.es!antodo	55	ierr = VecAXPY(st->w,-st->sigma,x);CHKERRQ(ierr);
778	dsic.upv.es!antodo	56	}
		57	ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
		58	if (st->sigma != 0.0) {
828	dsic.upv.es!antodo	59	ierr = VecAXPY(y,-st->sigma,st->w);CHKERRQ(ierr);
778	dsic.upv.es!antodo	60	}
		61	}
		62	PetscFunctionReturn(0);
		63	}
		64
		65	#undef __FUNCT__
813	dsic.upv.es!antodo	66	#define __FUNCT__ "STApplyTranspose_Fold"
		67	PetscErrorCode STApplyTranspose_Fold(ST st,Vec x,Vec y)
		68	{
		69	PetscErrorCode ierr;
		70	ST_FOLD ctx = (ST_FOLD ) st->data;
		71
		72	PetscFunctionBegin;
		73	if (st->B) {
		74	/* generalized eigenproblem: y = (A^T B^-T + sI)^2 x */
828	dsic.upv.es!antodo	75	ierr = STAssociatedKSPSolveTranspose(st,x,st->w);CHKERRQ(ierr);
		76	ierr = MatMult(st->A,st->w,ctx->w2);CHKERRQ(ierr);
813	dsic.upv.es!antodo	77	if (st->sigma != 0.0) {
828	dsic.upv.es!antodo	78	ierr = VecAXPY(ctx->w2,-st->sigma,x);CHKERRQ(ierr);
813	dsic.upv.es!antodo	79	}
828	dsic.upv.es!antodo	80	ierr = STAssociatedKSPSolveTranspose(st,ctx->w2,st->w);CHKERRQ(ierr);
		81	ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
		82	if (st->sigma != 0.0) {
		83	ierr = VecAXPY(y,-st->sigma,ctx->w2);CHKERRQ(ierr);
		84	}
813	dsic.upv.es!antodo	85	} else {
		86	/* standard eigenproblem: y = (A^T + sI)^2 x */
		87	ierr = MatMultTranspose(st->A,x,st->w);CHKERRQ(ierr);
		88	if (st->sigma != 0.0) {
828	dsic.upv.es!antodo	89	ierr = VecAXPY(st->w,-st->sigma,x);CHKERRQ(ierr);
813	dsic.upv.es!antodo	90	}
		91	ierr = MatMultTranspose(st->A,st->w,y);CHKERRQ(ierr);
		92	if (st->sigma != 0.0) {
828	dsic.upv.es!antodo	93	ierr = VecAXPY(y,-st->sigma,st->w);CHKERRQ(ierr);
813	dsic.upv.es!antodo	94	}
		95	}
		96	PetscFunctionReturn(0);
		97	}
		98
		99	#undef __FUNCT__
778	dsic.upv.es!antodo	100	#define __FUNCT__ "STBackTransform_Fold"
1779	antodo	101	PetscErrorCode STBackTransform_Fold(ST st,PetscInt n,PetscScalar eigr,PetscScalar eigi)
778	dsic.upv.es!antodo	102	{
1778	antodo	103	PetscInt j;
778	dsic.upv.es!antodo	104	PetscFunctionBegin;
1778	antodo	105	PetscValidScalarPointer(eigr,3);
		106	PetscValidScalarPointer(eigi,4);
		107	for (j=0;j<n;j++) {
778	dsic.upv.es!antodo	108	#if !defined(PETSC_USE_COMPLEX)
1778	antodo	109	if (eigi[j] == 0) {
778	dsic.upv.es!antodo	110	#endif
2089	jroman	111	eigr[j] = st->sigma + PetscSqrtScalar(eigr[j]);
778	dsic.upv.es!antodo	112	#if !defined(PETSC_USE_COMPLEX)
		113	} else {
1778	antodo	114	PetscScalar r,x,y;
		115	r = PetscSqrtScalar(eigr[j] * eigr[j] + eigi[j] * eigi[j]);
		116	x = PetscSqrtScalar((r + eigr[j]) / 2);
		117	y = PetscSqrtScalar((r - eigr[j]) / 2);
		118	if (eigi[j] < 0) y = - y;
2089	jroman	119	eigr[j] = st->sigma + x;
		120	eigi[j] = y;
778	dsic.upv.es!antodo	121	}
1778	antodo	122	#endif
778	dsic.upv.es!antodo	123	}
		124	PetscFunctionReturn(0);
		125	}
		126
		127	#undef __FUNCT__
		128	#define __FUNCT__ "STSetUp_Fold"
		129	PetscErrorCode STSetUp_Fold(ST st)
		130	{
		131	PetscErrorCode ierr;
		132	ST_FOLD ctx = (ST_FOLD ) st->data;
		133
		134	PetscFunctionBegin;
2087	jroman	135	/* if the user did not set the shift, use the target value */
		136	if (!st->sigma_set) st->sigma = st->defsigma;
		137
778	dsic.upv.es!antodo	138	if (st->B) {
		139	ierr = KSPSetOperators(st->ksp,st->B,st->B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
		140	ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
		141	if (ctx->w2) { ierr = VecDestroy(ctx->w2);CHKERRQ(ierr); }
		142	ierr = MatGetVecs(st->B,&ctx->w2,PETSC_NULL);CHKERRQ(ierr);
		143	}
		144	PetscFunctionReturn(0);
		145	}
		146
		147
		148	#undef __FUNCT__
		149	#define __FUNCT__ "STSetFromOptions_Fold"
		150	PetscErrorCode STSetFromOptions_Fold(ST st)
		151	{
		152	PetscErrorCode ierr;
2005	eromero	153	PC pc;
		154	const PCType pctype;
		155	const KSPType ksptype;
778	dsic.upv.es!antodo	156
		157	PetscFunctionBegin;
2005	eromero	158
		159	ierr = KSPGetPC(st->ksp,&pc);CHKERRQ(ierr);
		160	ierr = KSPGetType(st->ksp,&ksptype);CHKERRQ(ierr);
		161	ierr = PCGetType(pc,&pctype);CHKERRQ(ierr);
		162	if (!pctype && !ksptype) {
		163	if (st->shift_matrix == ST_MATMODE_SHELL) {
		164	/* in shell mode use GMRES with Jacobi as the default */
		165	ierr = KSPSetType(st->ksp,KSPGMRES);CHKERRQ(ierr);
		166	ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
		167	} else {
		168	/* use direct solver as default */
		169	ierr = KSPSetType(st->ksp,KSPPREONLY);CHKERRQ(ierr);
		170	ierr = PCSetType(pc,PCREDUNDANT);CHKERRQ(ierr);
		171	}
		172	}
		173
778	dsic.upv.es!antodo	174	PetscFunctionReturn(0);
		175	}
		176
		177	#undef __FUNCT__
		178	#define __FUNCT__ "STDestroy_Fold"
		179	PetscErrorCode STDestroy_Fold(ST st)
		180	{
		181	PetscErrorCode ierr;
		182	ST_FOLD ctx = (ST_FOLD ) st->data;
		183
		184	PetscFunctionBegin;
		185	if (ctx->w2) { ierr = VecDestroy(ctx->w2);CHKERRQ(ierr); }
		186	ierr = PetscFree(ctx);CHKERRQ(ierr);
		187	PetscFunctionReturn(0);
		188	}
		189
		190	EXTERN_C_BEGIN
		191	#undef __FUNCT__
		192	#define __FUNCT__ "STCreate_Fold"
		193	PetscErrorCode STCreate_Fold(ST st)
		194	{
		195	PetscErrorCode ierr;
		196	ST_FOLD *ctx;
		197
		198	PetscFunctionBegin;
		199
		200	ierr = PetscNew(ST_FOLD,&ctx); CHKERRQ(ierr);
		201	PetscLogObjectMemory(st,sizeof(ST_FOLD));
		202	st->data = (void *) ctx;
		203
1358	slepc	204	st->ops->apply = STApply_Fold;
		205	st->ops->getbilinearform = STGetBilinearForm_Default;
		206	st->ops->applytrans = STApplyTranspose_Fold;
		207	st->ops->backtr = STBackTransform_Fold;
		208	st->ops->setup = STSetUp_Fold;
2089	jroman	209	st->ops->view = STView_Default;
1358	slepc	210	st->ops->setfromoptions = STSetFromOptions_Fold;
		211	st->ops->destroy = STDestroy_Fold;
		212	st->checknullspace = 0;
778	dsic.upv.es!antodo	213
		214	PetscFunctionReturn(0);
		215	}
		216	EXTERN_C_END
		217

Subversion Repositories slepc-dev

(root)/trunk/src/st/impls/fold/fold.c @ 2216 - Rev 2116