WebSVN - slepc-dev - Blame - Rev 476 - /trunk/src/st/impls/cayley/cayley.c

Rev	Author	Line No.	Line
352	dsic.upv.es!antodo	1	/*
354	dsic.upv.es!jroman	2	Implements the Cayley spectral transform.
352	dsic.upv.es!antodo	3	*/
		4	#include "src/st/stimpl.h" /I "slepcst.h" I/
		5
		6	typedef struct {
		7	PetscScalar tau;
359	dsic.upv.es!antodo	8	PetscTruth tau_set;
352	dsic.upv.es!antodo	9	Vec w2;
		10	} ST_CAYLEY;
		11
		12	#undef __FUNCT__
		13	#define __FUNCT__ "STApply_Cayley"
476	dsic.upv.es!antodo	14	PetscErrorCode STApply_Cayley(ST st,Vec x,Vec y)
352	dsic.upv.es!antodo	15	{
476	dsic.upv.es!antodo	16	PetscErrorCode ierr;
		17	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
		18	PetscScalar tau = ctx->tau;
352	dsic.upv.es!antodo	19
354	dsic.upv.es!jroman	20	PetscFunctionBegin;
352	dsic.upv.es!antodo	21	if (st->shift_matrix == STMATMODE_INPLACE) { tau = tau + st->sigma; };
		22
		23	if (st->B) {
		24	/* generalized eigenproblem: y = (A - sB)^-1 (A + tB)x */
		25	ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
		26	ierr = MatMult(st->B,x,ctx->w2);CHKERRQ(ierr);
		27	ierr = VecAXPY(&tau,ctx->w2,st->w);CHKERRQ(ierr);
		28	ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
		29	}
		30	else {
		31	/* standard eigenproblem: y = (A - sI)^-1 (A + tI)x */
		32	ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
		33	ierr = VecAXPY(&tau,x,st->w);CHKERRQ(ierr);
		34	ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
		35	}
		36	PetscFunctionReturn(0);
		37	}
		38
		39	#undef __FUNCT__
		40	#define __FUNCT__ "STApplyNoB_Cayley"
476	dsic.upv.es!antodo	41	PetscErrorCode STApplyNoB_Cayley(ST st,Vec x,Vec y)
352	dsic.upv.es!antodo	42	{
476	dsic.upv.es!antodo	43	PetscErrorCode ierr;
352	dsic.upv.es!antodo	44
		45	PetscFunctionBegin;
		46	ierr = STAssociatedKSPSolve(st,x,y);CHKERRQ(ierr);
		47	PetscFunctionReturn(0);
		48	}
		49
		50	#undef __FUNCT__
		51	#define __FUNCT__ "STApplyB_Cayley"
476	dsic.upv.es!antodo	52	PetscErrorCode STApplyB_Cayley(ST st,Vec x,Vec y)
352	dsic.upv.es!antodo	53	{
476	dsic.upv.es!antodo	54	PetscErrorCode ierr;
		55	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
		56	PetscScalar tau = ctx->tau;
352	dsic.upv.es!antodo	57
354	dsic.upv.es!jroman	58	PetscFunctionBegin;
352	dsic.upv.es!antodo	59	if (st->shift_matrix == STMATMODE_INPLACE) { tau = tau + st->sigma; };
		60
		61	if (st->B) {
		62	/* generalized eigenproblem: y = (A + tB)x */
354	dsic.upv.es!jroman	63	ierr = MatMult(st->A,x,y);CHKERRQ(ierr);
352	dsic.upv.es!antodo	64	ierr = MatMult(st->B,x,ctx->w2);CHKERRQ(ierr);
354	dsic.upv.es!jroman	65	ierr = VecAXPY(&tau,ctx->w2,y);CHKERRQ(ierr);
352	dsic.upv.es!antodo	66	}
		67	else {
		68	/* standard eigenproblem: y = (A + tI)x */
354	dsic.upv.es!jroman	69	ierr = MatMult(st->A,x,y);CHKERRQ(ierr);
		70	ierr = VecAXPY(&tau,x,y);CHKERRQ(ierr);
352	dsic.upv.es!antodo	71	}
		72	PetscFunctionReturn(0);
		73	}
		74
		75	#undef __FUNCT__
		76	#define __FUNCT__ "STBackTransform_Cayley"
476	dsic.upv.es!antodo	77	PetscErrorCode STBackTransform_Cayley(ST st,PetscScalar eigr,PetscScalar eigi)
352	dsic.upv.es!antodo	78	{
		79	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
		80	#ifndef PETSC_USE_COMPLEX
		81	PetscScalar t,i,r;
		82	PetscFunctionBegin;
		83	PetscValidPointer(eigr,2);
		84	PetscValidPointer(eigi,3);
354	dsic.upv.es!jroman	85	if (eigi == 0.0) eigr = (ctx->tau + eigr st->sigma) / (*eigr - 1.0);
352	dsic.upv.es!antodo	86	else {
		87	r = *eigr;
		88	i = *eigi;
		89	r = st->sigma * (r * r + i * i - r) + ctx->tau * (r - 1);
		90	i = - st->sigma * i - ctx->tau * i;
		91	t = i * i + r * (r - 2.0) + 1.0;
		92	*eigr = r / t;
		93	*eigi = i / t;
		94	}
		95	#else
		96	PetscFunctionBegin;
		97	PetscValidPointer(eigr,2);
		98	eigr = (ctx->tau + eigr * st->sigma) / (*eigr - 1.0);
		99	#endif
		100	PetscFunctionReturn(0);
		101	}
		102
		103	#undef __FUNCT__
		104	#define __FUNCT__ "STPost_Cayley"
476	dsic.upv.es!antodo	105	PetscErrorCode STPost_Cayley(ST st)
352	dsic.upv.es!antodo	106	{
476	dsic.upv.es!antodo	107	PetscErrorCode ierr;
		108	PetscScalar alpha;
352	dsic.upv.es!antodo	109
		110	PetscFunctionBegin;
		111	if (st->shift_matrix == STMATMODE_INPLACE) {
		112	alpha = st->sigma;
		113	if( st->B ) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
		114	else { ierr = MatShift(&alpha, st->A); CHKERRQ(ierr); }
		115	st->setupcalled = 0;
		116	}
		117	PetscFunctionReturn(0);
		118	}
		119
		120	#undef __FUNCT__
		121	#define __FUNCT__ "STSetUp_Cayley"
476	dsic.upv.es!antodo	122	PetscErrorCode STSetUp_Cayley(ST st)
352	dsic.upv.es!antodo	123	{
476	dsic.upv.es!antodo	124	PetscErrorCode ierr;
		125	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
		126	PetscScalar alpha;
352	dsic.upv.es!antodo	127
		128	PetscFunctionBegin;
		129
359	dsic.upv.es!antodo	130	if (st->mat) { ierr = MatDestroy(st->mat);CHKERRQ(ierr); }
		131	if (!ctx->tau_set) { ctx->tau = st->sigma; }
		132	if (ctx->tau == 0.0 && st->sigma == 0.0) {
		133	SETERRQ(1,"Values of shift and antishift cannot be zero simultaneously");
		134	}
		135
352	dsic.upv.es!antodo	136	switch (st->shift_matrix) {
		137	case STMATMODE_INPLACE:
359	dsic.upv.es!antodo	138	st->mat = PETSC_NULL;
352	dsic.upv.es!antodo	139	if (st->sigma != 0.0) {
		140	alpha = -st->sigma;
		141	if (st->B) {
		142	ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr);
		143	} else {
		144	ierr = MatShift(&alpha,st->A);CHKERRQ(ierr);
		145	}
		146	}
		147	/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
		148	* improve performance when solving a number of related eigenproblems */
		149	ierr = KSPSetOperators(st->ksp,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
		150	break;
		151	case STMATMODE_SHELL:
		152	ierr = STMatShellCreate(st,&st->mat);CHKERRQ(ierr);
		153	ierr = KSPSetOperators(st->ksp,st->mat,st->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
		154	break;
		155	default:
		156	ierr = MatDuplicate(st->A,MAT_COPY_VALUES,&st->mat);CHKERRQ(ierr);
		157	if (st->sigma != 0.0) {
		158	alpha = -st->sigma;
		159	if (st->B) {
		160	ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr);
		161	} else {
		162	ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr);
		163	}
		164	}
		165	/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
		166	* improve performance when solving a number of related eigenproblems */
		167	ierr = KSPSetOperators(st->ksp,st->mat,st->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
		168	}
		169	if (st->B) {
		170	if (ctx->w2) { ierr = VecDestroy(ctx->w2);CHKERRQ(ierr); }
		171	ierr = MatGetVecs(st->B,&ctx->w2,PETSC_NULL);CHKERRQ(ierr);
359	dsic.upv.es!antodo	172	}
352	dsic.upv.es!antodo	173	ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
		174	PetscFunctionReturn(0);
		175	}
		176
		177	#undef __FUNCT__
		178	#define __FUNCT__ "STSetShift_Cayley"
476	dsic.upv.es!antodo	179	PetscErrorCode STSetShift_Cayley(ST st,PetscScalar newshift)
352	dsic.upv.es!antodo	180	{
476	dsic.upv.es!antodo	181	PetscErrorCode ierr;
		182	PetscScalar alpha;
		183	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
352	dsic.upv.es!antodo	184
		185	PetscFunctionBegin;
		186
		187	/* Nothing to be done if STSetUp has not been called yet */
		188	if (!st->setupcalled) PetscFunctionReturn(0);
		189
		190	switch (st->shift_matrix) {
		191	case STMATMODE_INPLACE:
		192	/* Undo previous operations */
		193	if (st->sigma != 0.0) {
		194	alpha = st->sigma;
		195	if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
		196	else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
		197	}
		198	/* Apply new shift */
		199	if (newshift != 0.0) {
		200	alpha = -newshift;
		201	if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
		202	else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
		203	}
		204	ierr = KSPSetOperators(st->ksp,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
		205	break;
		206	case STMATMODE_SHELL:
354	dsic.upv.es!jroman	207	ierr = KSPSetOperators(st->ksp,st->mat,st->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
352	dsic.upv.es!antodo	208	break;
		209	default:
354	dsic.upv.es!jroman	210	ierr = MatCopy(st->A, st->mat,SUBSET_NONZERO_PATTERN); CHKERRQ(ierr);
352	dsic.upv.es!antodo	211	if (newshift != 0.0) {
		212	alpha = -newshift;
		213	if (st->B) { ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr); }
		214	else { ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr); }
		215	}
		216	/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
		217	* improve performance when solving a number of related eigenproblems */
		218	ierr = KSPSetOperators(st->ksp,st->mat,st->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
		219	}
		220	st->sigma = newshift;
359	dsic.upv.es!antodo	221	if (!ctx->tau_set) { ctx->tau = newshift; }
352	dsic.upv.es!antodo	222	ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
		223	PetscFunctionReturn(0);
		224	}
		225
		226	#undef __FUNCT__
359	dsic.upv.es!antodo	227	#define __FUNCT__ "STSetFromOptions_Cayley"
476	dsic.upv.es!antodo	228	PetscErrorCode STSetFromOptions_Cayley(ST st)
		229	{
		230	PetscErrorCode ierr;
		231	PetscScalar tau;
		232	PetscTruth flg;
		233	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
359	dsic.upv.es!antodo	234
		235	PetscFunctionBegin;
		236	ierr = PetscOptionsHead("ST Cayley Options");CHKERRQ(ierr);
		237	ierr = PetscOptionsScalar("-st_antishift","Value of the antishift","STSetAntishift",ctx->tau,&tau,&flg); CHKERRQ(ierr);
		238	if (flg) {
		239	ierr = STCayleySetAntishift(st,tau);CHKERRQ(ierr);
		240	}
		241	ierr = PetscOptionsTail();CHKERRQ(ierr);
		242	PetscFunctionReturn(0);
		243	}
		244
		245	EXTERN_C_BEGIN
		246	#undef __FUNCT__
		247	#define __FUNCT__ "STCayleySetAntishift_Cayley"
476	dsic.upv.es!antodo	248	PetscErrorCode STCayleySetAntishift_Cayley(ST st,PetscScalar newshift)
359	dsic.upv.es!antodo	249	{
		250	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
		251
		252	PetscFunctionBegin;
		253	ctx->tau = newshift;
		254	ctx->tau_set = PETSC_TRUE;
		255	PetscFunctionReturn(0);
		256	}
		257	EXTERN_C_END
		258
		259	#undef __FUNCT__
352	dsic.upv.es!antodo	260	#define __FUNCT__ "STCayleySetAntishift"
476	dsic.upv.es!antodo	261	PetscErrorCode STCayleySetAntishift(ST st,PetscScalar newshift)
352	dsic.upv.es!antodo	262	{
476	dsic.upv.es!antodo	263	PetscErrorCode ierr, (*f)(ST,PetscScalar);
359	dsic.upv.es!antodo	264
		265	PetscFunctionBegin;
		266	PetscValidHeaderSpecific(st,ST_COOKIE,1);
		267	ierr = PetscObjectQueryFunction((PetscObject)st,"STCayleySetAntishift_C",(void (**)(void))&f);CHKERRQ(ierr);
		268	if (f) {
		269	ierr = (*f)(st,newshift);CHKERRQ(ierr);
		270	}
		271	PetscFunctionReturn(0);
		272	}
		273
		274	#undef __FUNCT__
		275	#define __FUNCT__ "STView_Cayley"
476	dsic.upv.es!antodo	276	PetscErrorCode STView_Cayley(ST st,PetscViewer viewer)
359	dsic.upv.es!antodo	277	{
476	dsic.upv.es!antodo	278	PetscErrorCode ierr;
		279	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
352	dsic.upv.es!antodo	280
		281	PetscFunctionBegin;
359	dsic.upv.es!antodo	282	#if !defined(PETSC_USE_COMPLEX)
		283	ierr = PetscViewerASCIIPrintf(viewer," antishift: %g\n",ctx->tau);CHKERRQ(ierr);
		284	#else
		285	ierr = PetscViewerASCIIPrintf(viewer," antishift: %g+%g i\n",PetscRealPart(ctx->tau),PetscImaginaryPart(ctx->tau));CHKERRQ(ierr);
		286	#endif
438	dsic.upv.es!antodo	287	ierr = STView_Default(st,viewer);CHKERRQ(ierr);
352	dsic.upv.es!antodo	288	PetscFunctionReturn(0);
		289	}
		290
		291	#undef __FUNCT__
		292	#define __FUNCT__ "STDestroy_Cayley"
476	dsic.upv.es!antodo	293	PetscErrorCode STDestroy_Cayley(ST st)
352	dsic.upv.es!antodo	294	{
476	dsic.upv.es!antodo	295	PetscErrorCode ierr;
		296	ST_CAYLEY ctx = (ST_CAYLEY ) st->data;
352	dsic.upv.es!antodo	297
		298	PetscFunctionBegin;
		299	if (ctx->w2) { ierr = VecDestroy(ctx->w2);CHKERRQ(ierr); }
		300	ierr = PetscFree(ctx);CHKERRQ(ierr);
		301	PetscFunctionReturn(0);
		302	}
		303
		304	EXTERN_C_BEGIN
		305	#undef __FUNCT__
		306	#define __FUNCT__ "STCreate_Cayley"
476	dsic.upv.es!antodo	307	PetscErrorCode STCreate_Cayley(ST st)
352	dsic.upv.es!antodo	308	{
476	dsic.upv.es!antodo	309	PetscErrorCode ierr;
		310	ST_CAYLEY *ctx;
352	dsic.upv.es!antodo	311
		312	PetscFunctionBegin;
		313	ierr = PetscNew(ST_CAYLEY,&ctx); CHKERRQ(ierr);
		314	PetscMemzero(ctx,sizeof(ST_CAYLEY));
		315	PetscLogObjectMemory(st,sizeof(ST_CAYLEY));
		316	st->data = (void *) ctx;
		317
		318	st->ops->apply = STApply_Cayley;
		319	st->ops->applyB = STApplyB_Cayley;
		320	st->ops->applynoB = STApplyNoB_Cayley;
		321	st->ops->postsolve = STPost_Cayley;
		322	st->ops->backtr = STBackTransform_Cayley;
		323	st->ops->setfromoptions = STSetFromOptions_Cayley;
		324	st->ops->setup = STSetUp_Cayley;
		325	st->ops->setshift = STSetShift_Cayley;
		326	st->ops->destroy = STDestroy_Cayley;
359	dsic.upv.es!antodo	327	st->ops->view = STView_Cayley;
352	dsic.upv.es!antodo	328
		329	st->checknullspace = STCheckNullSpace_Default;
		330
		331	ctx->tau = 0.0;
359	dsic.upv.es!antodo	332	ctx->tau_set = PETSC_FALSE;
352	dsic.upv.es!antodo	333
359	dsic.upv.es!antodo	334	ierr = PetscObjectComposeFunctionDynamic((PetscObject)st,"STCayleySetAntishift_C","STCayleySetAntishift_Cayley",
		335	STCayleySetAntishift_Cayley);CHKERRQ(ierr);
		336
352	dsic.upv.es!antodo	337	PetscFunctionReturn(0);
		338	}
		339	EXTERN_C_END
		340

Subversion Repositories slepc-dev

(root)/trunk/src/st/impls/cayley/cayley.c @ 1672 - Rev 476