Subversion Repositories slepc-dev

/*

     SVD routines for setting up the solver.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*/

#include "private/svdimpl.h"      /*I "slepcsvd.h" I*/

#undef __FUNCT__  

#define __FUNCT__ "SVDSetOperator"

/*@

   SVDSetOperator - Set the matrix associated with the singular value problem.

   Collective on SVD and Mat

   Input Parameters:

+  svd - the singular value solver context

-  A  - the matrix associated with the singular value problem

   Level: beginner

.seealso: SVDSolve(), SVDGetOperator()

@*/

PetscErrorCode SVDSetOperator(SVD svd,Mat mat)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  PetscValidHeaderSpecific(mat,MAT_COOKIE,2);

  PetscCheckSameComm(svd,1,mat,2);

  ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr);

  if (svd->OP) {

    ierr = MatDestroy(svd->OP);CHKERRQ(ierr);

  }

  svd->OP = mat;

  if (svd->vec_initial) {

    ierr = VecDestroy(svd->vec_initial);CHKERRQ(ierr);

    svd->vec_initial = PETSC_NULL;

  }

  svd->setupcalled = 0;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "SVDGetOperator"

/*@

   SVDGetOperator - Get the matrix associated with the singular value problem.

   Not collective, though parallel Mats are returned if the SVD is parallel

   Input Parameter:

.  svd - the singular value solver context

   Output Parameters:

.  A    - the matrix associated with the singular value problem

   Level: advanced

.seealso: SVDSolve(), SVDSetOperator()

@*/

PetscErrorCode SVDGetOperator(SVD svd,Mat *A)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  PetscValidPointer(A,2);

  *A = svd->OP;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "SVDSetInitialVector"

/*@

   SVDSetInitialVector - Sets the initial vector from which the

   singular value solver starts to iterate.

   Collective on SVD and Vec

   Input Parameters:

+  svd - the singular value solver context

-  vec - the vector

   Level: intermediate

.seealso: SVDGetInitialVector()

@*/

PetscErrorCode SVDSetInitialVector(SVD svd,Vec vec)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  PetscValidHeaderSpecific(vec,VEC_COOKIE,2);

  PetscCheckSameComm(svd,1,vec,2);

  ierr = PetscObjectReference((PetscObject)vec);CHKERRQ(ierr);

  if (svd->vec_initial) {

    ierr = VecDestroy(svd->vec_initial); CHKERRQ(ierr);

  }

  svd->vec_initial = vec;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "SVDGetInitialVector"

/*@

   SVDGetInitialVector - Gets the initial vector associated with the

   singular value solver; if the vector was not set it will return a 0

   pointer or a vector randomly generated by SVDSetUp().

   Not collective, but vector is shared by all processors that share the SVD

   Input Parameter:

.  svd - the singular value solver context

   Output Parameter:

.  vec - the vector

   Level: intermediate

.seealso: SVDSetInitialVector()

@*/

PetscErrorCode SVDGetInitialVector(SVD svd,Vec *vec)

{

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  PetscValidPointer(vec,2);

  *vec = svd->vec_initial;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "SVDSetUp"

/*@

   SVDSetUp - Sets up all the internal data structures necessary for the

   execution of the singular value solver.

   Collective on SVD

   Input Parameter:

.  svd   - singular value solver context

   Level: advanced

   Notes:

   This function need not be called explicitly in most cases, since SVDSolve()

   calls it. It can be useful when one wants to measure the set-up time

   separately from the solve time.

.seealso: SVDCreate(), SVDSolve(), SVDDestroy()

@*/

PetscErrorCode SVDSetUp(SVD svd)

{

  PetscErrorCode ierr;

  PetscTruth     flg;

  PetscInt       i,M,N,nloc;

  PetscScalar    *pV;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  if (svd->setupcalled) PetscFunctionReturn(0);

  ierr = PetscLogEventBegin(SVD_SetUp,svd,0,0,0);CHKERRQ(ierr);

  /* Set default solver type */

  if (!((PetscObject)svd)->type_name) {

    ierr = SVDSetType(svd,SVDCROSS);CHKERRQ(ierr);

  }

  /* check matrix */

  if (!svd->OP)

    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSetOperator must be called first");

  /* determine how to build the transpose */

  if (svd->transmode == PETSC_DECIDE) {

    ierr = MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg);CHKERRQ(ierr);    

    if (flg) svd->transmode = SVD_TRANSPOSE_EXPLICIT;

    else svd->transmode = SVD_TRANSPOSE_IMPLICIT;

  }

  /* build transpose matrix */

  if (svd->A) { ierr = MatDestroy(svd->A);CHKERRQ(ierr); }

  if (svd->AT) { ierr = MatDestroy(svd->AT);CHKERRQ(ierr); }

  ierr = MatGetSize(svd->OP,&M,&N);CHKERRQ(ierr);

  ierr = PetscObjectReference((PetscObject)svd->OP);CHKERRQ(ierr);

  switch (svd->transmode) {

    case SVD_TRANSPOSE_EXPLICIT:

      ierr = MatHasOperation(svd->OP,MATOP_TRANSPOSE,&flg);CHKERRQ(ierr);

      if (!flg) SETERRQ(1,"Matrix has not defined the MatTranpose operation");

      if (M>=N) {

        svd->A = svd->OP;

        ierr = MatTranspose(svd->OP, MAT_INITIAL_MATRIX,&svd->AT);CHKERRQ(ierr);

      } else {

        ierr = MatTranspose(svd->OP, MAT_INITIAL_MATRIX,&svd->A);CHKERRQ(ierr);

        svd->AT = svd->OP;

      }

      break;

    case SVD_TRANSPOSE_IMPLICIT:

      ierr = MatHasOperation(svd->OP,MATOP_MULT_TRANSPOSE,&flg);CHKERRQ(ierr);

      if (!flg) SETERRQ(1,"Matrix has not defined the MatMultTranpose operation");

      if (M>=N) {

        svd->A = svd->OP;

        svd->AT = PETSC_NULL;    

      } else {

        svd->A = PETSC_NULL;

        svd->AT = svd->OP;

      }

      break;

    default:

      SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid transpose mode");

  }

  /* set initial vector */

  if (!svd->vec_initial) {

    ierr = SVDMatGetVecs(svd,&svd->vec_initial,PETSC_NULL);CHKERRQ(ierr);

    ierr = SlepcVecSetRandom(svd->vec_initial);CHKERRQ(ierr);

  }

  /* call specific solver setup */

  ierr = (*svd->ops->setup)(svd);CHKERRQ(ierr);

  if (svd->ncv > M || svd->ncv > N)

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"ncv bigger than matrix dimensions");

  if (svd->nsv > svd->ncv)

    SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"nsv bigger than ncv");

  if (svd->ncv != svd->n) {

    /* free memory for previous solution  */

    if (svd->n) {

      ierr = PetscFree(svd->sigma);CHKERRQ(ierr);

      ierr = PetscFree(svd->perm);CHKERRQ(ierr);

      ierr = PetscFree(svd->errest);CHKERRQ(ierr);

      ierr = VecGetArray(svd->V[0],&pV);CHKERRQ(ierr);

      for (i=0;i<svd->n;i++) {

        ierr = VecDestroy(svd->V[i]);CHKERRQ(ierr);

      }

      ierr = PetscFree(pV);CHKERRQ(ierr);

      ierr = PetscFree(svd->V);CHKERRQ(ierr);

    }

    /* allocate memory for next solution */

    ierr = PetscMalloc(svd->ncv*sizeof(PetscReal),&svd->sigma);CHKERRQ(ierr);

    ierr = PetscMalloc(svd->ncv*sizeof(PetscReal),&svd->perm);CHKERRQ(ierr);

    ierr = PetscMalloc(svd->ncv*sizeof(PetscReal),&svd->errest);CHKERRQ(ierr);

    ierr = PetscMalloc(svd->ncv*sizeof(Vec),&svd->V);CHKERRQ(ierr);

    ierr = VecGetLocalSize(svd->vec_initial,&nloc);CHKERRQ(ierr);

    ierr = PetscMalloc(svd->ncv*nloc*sizeof(PetscScalar),&pV);CHKERRQ(ierr);

    for (i=0;i<svd->ncv;i++) {

      ierr = VecCreateMPIWithArray(((PetscObject)svd)->comm,nloc,PETSC_DECIDE,pV+i*nloc,&svd->V[i]);CHKERRQ(ierr);

    }

    svd->n = svd->ncv;

  }

  ierr = PetscLogEventEnd(SVD_SetUp,svd,0,0,0);CHKERRQ(ierr);

  svd->setupcalled = 1;

  PetscFunctionReturn(0);

}