WebSVN - slepc-dev - Diff - Rev 470 and 476 - /trunk/src/eps/impls/subspace/subspace.c

Rev 470	Rev 476
Line 2...	Line 2...
#include "src/eps/epsimpl.h" /I "slepceps.h" I/	#include "src/eps/epsimpl.h" /I "slepceps.h" I/
#include "slepcblaslapack.h"	#include "slepcblaslapack.h"

#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSSetUp_SUBSPACE"	#define __FUNCT__ "EPSSetUp_SUBSPACE"
static int EPSSetUp_SUBSPACE(EPS eps)	PetscErrorCode EPSSetUp_SUBSPACE(EPS eps)
{	{
int ierr, N;	PetscErrorCode ierr;
	int N;

PetscFunctionBegin;	PetscFunctionBegin;
ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);	ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);
if (eps->ncv) {	if (eps->ncv) {
if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");	if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");
Line 24...	Line 25...
PetscFunctionReturn(0);	PetscFunctionReturn(0);
}	}

#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSdcond"	#define __FUNCT__ "EPSdcond"
static int EPSdcond(PetscScalar* H,int n, PetscReal* cond)	static PetscErrorCode EPSdcond(PetscScalar* H,int n, PetscReal* cond)
{	{
int ierr,*ipiv,lwork;	PetscErrorCode ierr;
PetscScalar *work;	int *ipiv,lwork,info;
PetscReal hn,hin,*rwork;	PetscScalar *work;
	PetscReal hn,hin,*rwork;

PetscFunctionBegin;	PetscFunctionBegin;
ierr = PetscMalloc(sizeof(int)*n,&ipiv);CHKERRQ(ierr);	ierr = PetscMalloc(sizeof(int)*n,&ipiv);CHKERRQ(ierr);
lwork = n*n;	lwork = n*n;
ierr = PetscMalloc(sizeof(PetscScalar)*lwork,&work);CHKERRQ(ierr);	ierr = PetscMalloc(sizeof(PetscScalar)*lwork,&work);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscReal)*n,&rwork);CHKERRQ(ierr);	ierr = PetscMalloc(sizeof(PetscReal)*n,&rwork);CHKERRQ(ierr);
hn = LAlanhs_("I",&n,H,&n,rwork,1);	hn = LAlanhs_("I",&n,H,&n,rwork,1);
LAgetrf_(&n,&n,H,&n,ipiv,&ierr);	LAgetrf_(&n,&n,H,&n,ipiv,&info);
if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xGETRF");	if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGETRF %d",info);
LAgetri_(&n,H,&n,ipiv,work,&lwork,&ierr);	LAgetri_(&n,H,&n,ipiv,work,&lwork,&info);
if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xGETRI");	if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGETRI %d",info);
hin = LAlange_("I",&n,&n,H,&n,rwork,1);	hin = LAlange_("I",&n,&n,H,&n,rwork,1);
cond = hn hin;	cond = hn hin;
ierr = PetscFree(ipiv);CHKERRQ(ierr);	ierr = PetscFree(ipiv);CHKERRQ(ierr);
ierr = PetscFree(work);CHKERRQ(ierr);	ierr = PetscFree(work);CHKERRQ(ierr);
ierr = PetscFree(rwork);CHKERRQ(ierr);	ierr = PetscFree(rwork);CHKERRQ(ierr);
PetscFunctionReturn(0);	PetscFunctionReturn(0);
}	}

#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSdgroup"	#define __FUNCT__ "EPSdgroup"
static int EPSdgroup(int l,int m,PetscScalar wr,PetscScalar wi,PetscReal *rsd,	static PetscErrorCode EPSdgroup(int l,int m,PetscScalar wr,PetscScalar wi,PetscReal *rsd,
PetscReal grptol,int ngrp,PetscReal ctr,PetscReal ae,PetscReal arsd)	PetscReal grptol,int ngrp,PetscReal ctr,PetscReal ae,PetscReal arsd)
{	{
int i;	int i;
PetscReal rmod,rmod1;	PetscReal rmod,rmod1;

Line 99...	Line 101...
PetscFunctionReturn(0);	PetscFunctionReturn(0);
}	}

#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSSchurResidualNorms"	#define __FUNCT__ "EPSSchurResidualNorms"
static int EPSSchurResidualNorms(EPS eps,Vec V,Vec AV,PetscScalar T,int l,int m,int ldt,PetscReal rsd)	static PetscErrorCode EPSSchurResidualNorms(EPS eps,Vec V,Vec AV,PetscScalar T,int l,int m,int ldt,PetscReal rsd)
{	{
int ierr,i;	PetscErrorCode ierr;
PetscScalar zero = 0.0,minus = -1.0;	int i;
	PetscScalar zero = 0.0,minus = -1.0;
#if defined(PETSC_USE_COMPLEX)	#if defined(PETSC_USE_COMPLEX)
PetscScalar t;	PetscScalar t;
#endif	#endif

PetscFunctionBegin;	PetscFunctionBegin;
for (i=l;i<m;i++) {	for (i=l;i<m;i++) {
ierr = VecSet(&zero,eps->work[0]);CHKERRQ(ierr);	ierr = VecSet(&zero,eps->work[0]);CHKERRQ(ierr);
Line 136...	Line 139...
PetscFunctionReturn(0);	PetscFunctionReturn(0);
}	}

#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSSolve_SUBSPACE"	#define __FUNCT__ "EPSSolve_SUBSPACE"
static int EPSSolve_SUBSPACE(EPS eps)	PetscErrorCode EPSSolve_SUBSPACE(EPS eps)
{	{
int ierr,i,j,ilo,lwork,ngrp,nogrp,itrsd,itrsdold,	PetscErrorCode ierr;
nxtsrr,idsrr,*iwork,idort,nxtort,ncv = eps->ncv;	int i,j,ilo,lwork,info,ngrp,nogrp,itrsd,itrsdold,
PetscScalar T,U,tau,work,t;	nxtsrr,idsrr,*iwork,idort,nxtort,ncv = eps->ncv;
PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsdold,norm,tcond;	PetscScalar T,U,tau,work,t;
	PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsdold,norm,tcond;
/* Parameters */	/* Parameters */
int init = 5;	int init = 5;
PetscReal stpfac = 1.5,	PetscReal stpfac = 1.5,
alpha = 1.0,	alpha = 1.0,
beta = 1.1,	beta = 1.1,
grptol = 1e-8,	grptol = 1e-8,
cnvtol = 1e-6;	cnvtol = 1e-6;
int orttol = 2;	int orttol = 2;

PetscFunctionBegin;	PetscFunctionBegin;
eps->its = 0;	eps->its = 0;
eps->nconv = 0;	eps->nconv = 0;
ierr = PetscMalloc(sizeof(PetscScalar)ncvncv,&T);CHKERRQ(ierr);	ierr = PetscMalloc(sizeof(PetscScalar)ncvncv,&T);CHKERRQ(ierr);
Line 189...	Line 193...
ierr = VecMDot(ncv,eps->AV[i],eps->V,T+i*ncv);CHKERRQ(ierr);	ierr = VecMDot(ncv,eps->AV[i],eps->V,T+i*ncv);CHKERRQ(ierr);
}	}

/* [U,H] = hess(T) */	/* [U,H] = hess(T) */
ilo = eps->nconv + 1;	ilo = eps->nconv + 1;
LAgehrd_(&ncv,&ilo,&ncv,T,&ncv,tau,work,&lwork,&ierr);	LAgehrd_(&ncv,&ilo,&ncv,T,&ncv,tau,work,&lwork,&info);
if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xGEHRD");	if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGEHRD %d",info);
for (j=0;j<ncv-1;j++) {	for (j=0;j<ncv-1;j++) {
for (i=j+2;i<ncv;i++) {	for (i=j+2;i<ncv;i++) {
U[i+jncv] = T[i+jncv];	U[i+jncv] = T[i+jncv];
T[i+j*ncv] = 0.0;	T[i+j*ncv] = 0.0;
}	}
}	}
LAorghr_(&ncv,&ilo,&ncv,U,&ncv,tau,work,&lwork,&ierr);	LAorghr_(&ncv,&ilo,&ncv,U,&ncv,tau,work,&lwork,&info);
if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xORGHR");	if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xORGHR %d",info);

/* [T,wr,wi,U] = laqr3(H,U) */	/* [T,wr,wi,U] = laqr3(H,U) */
ierr = EPSDenseSchur(T,U,eps->eigr,eps->eigi,eps->nconv,ncv);CHKERRQ(ierr);	ierr = EPSDenseSchur(T,U,eps->eigr,eps->eigi,eps->nconv,ncv);CHKERRQ(ierr);

/* AV(:,idx) = AVU(:,idx) /	/* AV(:,idx) = AVU(:,idx) /
Line 309...	Line 313...
PetscFunctionReturn(0);	PetscFunctionReturn(0);
}	}

#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSComputeVectors_SUBSPACE"	#define __FUNCT__ "EPSComputeVectors_SUBSPACE"
int EPSComputeVectors_SUBSPACE(EPS eps)	PetscErrorCode EPSComputeVectors_SUBSPACE(EPS eps)
{	{
int ierr;	PetscErrorCode ierr;
PetscFunctionBegin;	PetscFunctionBegin;
ierr = EPSComputeVectors_Default(eps);CHKERRQ(ierr);	ierr = EPSComputeVectors_Default(eps);CHKERRQ(ierr);
PetscFunctionReturn(0);	PetscFunctionReturn(0);
}	}

EXTERN_C_BEGIN	EXTERN_C_BEGIN
#undef __FUNCT__	#undef __FUNCT__
#define __FUNCT__ "EPSCreate_SUBSPACE"	#define __FUNCT__ "EPSCreate_SUBSPACE"
int EPSCreate_SUBSPACE(EPS eps)	PetscErrorCode EPSCreate_SUBSPACE(EPS eps)
{	{
PetscFunctionBegin;	PetscFunctionBegin;
eps->ops->setup = EPSSetUp_SUBSPACE;	eps->ops->setup = EPSSetUp_SUBSPACE;
eps->ops->solve = EPSSolve_SUBSPACE;	eps->ops->solve = EPSSolve_SUBSPACE;
eps->ops->destroy = EPSDestroy_Default;	eps->ops->destroy = EPSDestroy_Default;

Line 2...

#include "src/eps/epsimpl.h"                /*I "slepceps.h" I*/

#include "src/eps/epsimpl.h"                /*I "slepceps.h" I*/

#include "slepcblaslapack.h"

#include "slepcblaslapack.h"

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSSetUp_SUBSPACE"

#define __FUNCT__ "EPSSetUp_SUBSPACE"

static int EPSSetUp_SUBSPACE(EPS eps)

PetscErrorCode EPSSetUp_SUBSPACE(EPS eps)

  int       ierr, N;

  PetscErrorCode ierr;

  int            N;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);

  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);

  if (eps->ncv) {

  if (eps->ncv) {

    if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");

    if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");

Line 24...

Line 25...

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSdcond"

#define __FUNCT__ "EPSdcond"

static int EPSdcond(PetscScalar* H,int n, PetscReal* cond)

static PetscErrorCode EPSdcond(PetscScalar* H,int n, PetscReal* cond)

  int         ierr,*ipiv,lwork;

  PetscErrorCode ierr;

  PetscScalar *work;

  int            *ipiv,lwork,info;

  PetscReal   hn,hin,*rwork;

  PetscScalar    *work;

  PetscReal      hn,hin,*rwork;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = PetscMalloc(sizeof(int)*n,&ipiv);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(int)*n,&ipiv);CHKERRQ(ierr);

  lwork = n*n;

  lwork = n*n;

  ierr = PetscMalloc(sizeof(PetscScalar)*lwork,&work);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscScalar)*lwork,&work);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscReal)*n,&rwork);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscReal)*n,&rwork);CHKERRQ(ierr);

  hn = LAlanhs_("I",&n,H,&n,rwork,1);

  hn = LAlanhs_("I",&n,H,&n,rwork,1);

  LAgetrf_(&n,&n,H,&n,ipiv,&ierr);

  LAgetrf_(&n,&n,H,&n,ipiv,&info);

  if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xGETRF");

  if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGETRF %d",info);

  LAgetri_(&n,H,&n,ipiv,work,&lwork,&ierr);

  LAgetri_(&n,H,&n,ipiv,work,&lwork,&info);

  if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xGETRI");

  if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGETRI %d",info);

  hin = LAlange_("I",&n,&n,H,&n,rwork,1);

  hin = LAlange_("I",&n,&n,H,&n,rwork,1);

  *cond = hn * hin;

  *cond = hn * hin;

  ierr = PetscFree(ipiv);CHKERRQ(ierr);

  ierr = PetscFree(ipiv);CHKERRQ(ierr);

  ierr = PetscFree(work);CHKERRQ(ierr);

  ierr = PetscFree(work);CHKERRQ(ierr);

  ierr = PetscFree(rwork);CHKERRQ(ierr);

  ierr = PetscFree(rwork);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSdgroup"

#define __FUNCT__ "EPSdgroup"

static int EPSdgroup(int l,int m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,

static PetscErrorCode EPSdgroup(int l,int m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,

  PetscReal grptol,int *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)

  PetscReal grptol,int *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)

  int       i;

  int       i;

  PetscReal rmod,rmod1;

  PetscReal rmod,rmod1;

Line 99...

Line 101...

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSSchurResidualNorms"

#define __FUNCT__ "EPSSchurResidualNorms"

static int EPSSchurResidualNorms(EPS eps,Vec *V,Vec *AV,PetscScalar *T,int l,int m,int ldt,PetscReal *rsd)

static PetscErrorCode EPSSchurResidualNorms(EPS eps,Vec *V,Vec *AV,PetscScalar *T,int l,int m,int ldt,PetscReal *rsd)

  int         ierr,i;

  PetscErrorCode ierr;

  PetscScalar zero = 0.0,minus = -1.0;

  int            i;

  PetscScalar    zero = 0.0,minus = -1.0;

#if defined(PETSC_USE_COMPLEX)

#if defined(PETSC_USE_COMPLEX)

  PetscScalar t;

  PetscScalar    t;

#endif

#endif

  PetscFunctionBegin;

  PetscFunctionBegin;

  for (i=l;i<m;i++) {

  for (i=l;i<m;i++) {

    ierr = VecSet(&zero,eps->work[0]);CHKERRQ(ierr);

    ierr = VecSet(&zero,eps->work[0]);CHKERRQ(ierr);

Line 136...

Line 139...

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSSolve_SUBSPACE"

#define __FUNCT__ "EPSSolve_SUBSPACE"

static int EPSSolve_SUBSPACE(EPS eps)

PetscErrorCode EPSSolve_SUBSPACE(EPS eps)

  int         ierr,i,j,ilo,lwork,ngrp,nogrp,*itrsd,*itrsdold,

  PetscErrorCode ierr;

              nxtsrr,idsrr,*iwork,idort,nxtort,ncv = eps->ncv;

  int            i,j,ilo,lwork,info,ngrp,nogrp,*itrsd,*itrsdold,

  PetscScalar *T,*U,*tau,*work,t;

                 nxtsrr,idsrr,*iwork,idort,nxtort,ncv = eps->ncv;

  PetscReal   arsd,oarsd,ctr,octr,ae,oae,*rsdold,norm,tcond;

  PetscScalar    *T,*U,*tau,*work,t;

  PetscReal      arsd,oarsd,ctr,octr,ae,oae,*rsdold,norm,tcond;

  /* Parameters */

  /* Parameters */

  int         init = 5;

  int            init = 5;

  PetscReal   stpfac = 1.5,

  PetscReal      stpfac = 1.5,

              alpha = 1.0,

                 alpha = 1.0,

              beta = 1.1,

                 beta = 1.1,

              grptol = 1e-8,

                 grptol = 1e-8,

              cnvtol = 1e-6;

                 cnvtol = 1e-6;

  int         orttol = 2;

  int            orttol = 2;

  PetscFunctionBegin;

  PetscFunctionBegin;

  eps->its = 0;

  eps->its = 0;

  eps->nconv = 0;

  eps->nconv = 0;

  ierr = PetscMalloc(sizeof(PetscScalar)*ncv*ncv,&T);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscScalar)*ncv*ncv,&T);CHKERRQ(ierr);

Line 189...

Line 193...

      ierr = VecMDot(ncv,eps->AV[i],eps->V,T+i*ncv);CHKERRQ(ierr);

      ierr = VecMDot(ncv,eps->AV[i],eps->V,T+i*ncv);CHKERRQ(ierr);

    /* [U,H] = hess(T) */

    /* [U,H] = hess(T) */

    ilo = eps->nconv + 1;

    ilo = eps->nconv + 1;

    LAgehrd_(&ncv,&ilo,&ncv,T,&ncv,tau,work,&lwork,&ierr);

    LAgehrd_(&ncv,&ilo,&ncv,T,&ncv,tau,work,&lwork,&info);

    if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xGEHRD");

    if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGEHRD %d",info);

    for (j=0;j<ncv-1;j++) {

    for (j=0;j<ncv-1;j++) {

      for (i=j+2;i<ncv;i++) {

      for (i=j+2;i<ncv;i++) {

        U[i+j*ncv] = T[i+j*ncv];

        U[i+j*ncv] = T[i+j*ncv];

        T[i+j*ncv] = 0.0;

        T[i+j*ncv] = 0.0;

    LAorghr_(&ncv,&ilo,&ncv,U,&ncv,tau,work,&lwork,&ierr);

    LAorghr_(&ncv,&ilo,&ncv,U,&ncv,tau,work,&lwork,&info);

    if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in Lapack xORGHR");

    if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xORGHR %d",info);

    /* [T,wr,wi,U] = laqr3(H,U) */

    /* [T,wr,wi,U] = laqr3(H,U) */

    ierr = EPSDenseSchur(T,U,eps->eigr,eps->eigi,eps->nconv,ncv);CHKERRQ(ierr);

    ierr = EPSDenseSchur(T,U,eps->eigr,eps->eigi,eps->nconv,ncv);CHKERRQ(ierr);

    /* AV(:,idx) = AV*U(:,idx) */

    /* AV(:,idx) = AV*U(:,idx) */

Line 309...

Line 313...

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSComputeVectors_SUBSPACE"

#define __FUNCT__ "EPSComputeVectors_SUBSPACE"

int EPSComputeVectors_SUBSPACE(EPS eps)

PetscErrorCode EPSComputeVectors_SUBSPACE(EPS eps)

  int ierr;

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = EPSComputeVectors_Default(eps);CHKERRQ(ierr);

  ierr = EPSComputeVectors_Default(eps);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

EXTERN_C_BEGIN

EXTERN_C_BEGIN

#undef __FUNCT__

#undef __FUNCT__

#define __FUNCT__ "EPSCreate_SUBSPACE"

#define __FUNCT__ "EPSCreate_SUBSPACE"

int EPSCreate_SUBSPACE(EPS eps)

PetscErrorCode EPSCreate_SUBSPACE(EPS eps)

  PetscFunctionBegin;

  PetscFunctionBegin;

  eps->ops->setup                = EPSSetUp_SUBSPACE;

  eps->ops->setup                = EPSSetUp_SUBSPACE;

  eps->ops->solve                = EPSSolve_SUBSPACE;

  eps->ops->solve                = EPSSolve_SUBSPACE;

  eps->ops->destroy              = EPSDestroy_Default;

  eps->ops->destroy              = EPSDestroy_Default;

Subversion Repositories slepc-dev

(root)/trunk/src/eps/impls/subspace/subspace.c @ 2329 - Rev 470 → 476