WebSVN - slepc-dev - Path Comparison - / Rev 1282 and / Rev 1283

Ignore whitespace Rev 1282 → Rev 1283

24,10 → 24,17

 #define SVDLAPACK      "lapack"
 #define SVDLANCZOS     "lanczos"
-typedef enum { SVD_TRANSPOSE_EXPLICIT, SVD_TRANSPOSE_MATMULT } SVDTransposeMode;
+typedef enum { SVD_TRANSPOSE_EXPLICIT, SVD_TRANSPOSE_IMPLICIT } SVDTransposeMode;
 typedef enum { SVD_LARGEST, SVD_SMALLEST } SVDWhich;
+typedef enum {/* converged */
+              SVD_CONVERGED_TOL                =  2,
+              /* diverged */
+              SVD_DIVERGED_ITS                 = -3,
+              SVD_DIVERGED_BREAKDOWN           = -4,
+              SVD_CONVERGED_ITERATING          =  0 } SVDConvergedReason;
 EXTERN PetscErrorCode SVDCreate(MPI_Comm,SVD*);
 EXTERN PetscErrorCode SVDSetType(SVD,SVDType);
 EXTERN PetscErrorCode SVDGetType(SVD,SVDType*);

46,6 → 53,7

EXTERN PetscErrorCode SVDSetFromOptions(SVD);

EXTERN PetscErrorCode SVDSetUp(SVD);

EXTERN PetscErrorCode SVDSolve(SVD);

EXTERN PetscErrorCode SVDGetIterationNumber(SVD,int*);

EXTERN PetscErrorCode SVDGetConverged(SVD,int*);

EXTERN PetscErrorCode SVDGetSingularTriplet(SVD,int,PetscReal*,Vec,Vec);

EXTERN PetscErrorCode SVDComputeResidualNorms(SVD,int,PetscReal*,PetscReal*);

53,7 → 61,7

 EXTERN PetscErrorCode SVDDestroy(SVD);
 EXTERN PetscErrorCode SVDInitializePackage(char*);
-typedef enum { SVDEIGENSOLVER_DIRECT, SVDEIGENSOLVER_TRANSPOSE,
+typedef enum { SVDEIGENSOLVER_ATA, SVDEIGENSOLVER_AAT,
                SVDEIGENSOLVER_CYCLIC } SVDEigensolverMode;
 EXTERN PetscErrorCode SVDEigensolverSetMode(SVD,SVDEigensolverMode);

/trunk/include/slepcblaslapack.h

74,7 → 74,8

 #define LAPACKgeevx_ SLEPC_BLASLAPACK(geevx,GEEVX)
 #define LAPACKggevx_ SLEPC_BLASLAPACK(ggevx,GGEVX)
 #define LAPACKgelqf_ SLEPC_BLASLAPACK(gelqf,GELQF)
-#define LAPACKgesdd  SLEPC_BLASLAPACK(gesdd,GESDD)
+#define LAPACKgesdd_ SLEPC_BLASLAPACK(gesdd,GESDD)
+#define LAPACKbdsqr_ SLEPC_BLASLAPACK(bdsqr,BDSQR)
 #if !defined(PETSC_USE_COMPLEX)
 #define LAPACKorghr_ SLEPC_BLASLAPACK(orghr,ORGHR)

104,6 → 105,7

EXTERN void LAPACKgelqf_(PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);

EXTERN void BLAStrsm_(const char*,const char*,const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt,PetscBLASInt,PetscBLASInt,PetscBLASInt);

EXTERN void LAPACKormlq_(const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);

EXTERN void LAPACKbdsqr_(const char*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscReal*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt);

#if !defined(PETSC_USE_COMPLEX)

EXTERN void LAPACKhseqr_(const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);

113,7 → 115,7

 EXTERN void      LAPACKggevx_(const char*,const char*,const char*,const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt,PetscBLASInt,PetscBLASInt);
 EXTERN void      LAPACKsyevr_(const char*,const char*,const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt,PetscBLASInt);
 EXTERN void      LAPACKsygvd_(PetscBLASInt*,const char*,const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);
-EXTERN void      LAPACKgesdd(const char*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt);
+EXTERN void      LAPACKgesdd_(const char*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt);
 #else
 EXTERN void      LAPACKhseqr_(const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);
 EXTERN void      LAPACKtrexc_(const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt);

122,7 → 124,7

 EXTERN void      LAPACKggevx_(const char*,const char*,const char*,const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*, PetscBLASInt*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscScalar*, PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt,PetscBLASInt,PetscBLASInt);
 EXTERN void      LAPACKsyevr_(const char *,const char*,const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscBLASInt*, PetscScalar*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt,PetscBLASInt);
 EXTERN void      LAPACKsygvd_(PetscBLASInt*,const char*,const char*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);
-EXTERN void      LAPACKgesdd(const char*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt);
+EXTERN void      LAPACKgesdd_(const char*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt);
 #endif

/trunk/src/svd/interface/svdbasic.c

112,12 → 112,21

       case SVD_TRANSPOSE_EXPLICIT:
         ierr = PetscViewerASCIIPrintf(viewer,"  transpose mode: explicit\n");CHKERRQ(ierr);
         break;
-      case SVD_TRANSPOSE_MATMULT:
-        ierr = PetscViewerASCIIPrintf(viewer,"  transpose mode: matmult\n");CHKERRQ(ierr);
+      case SVD_TRANSPOSE_IMPLICIT:
+        ierr = PetscViewerASCIIPrintf(viewer,"  transpose mode: implicit\n");CHKERRQ(ierr);
         break;
       default:
         ierr = PetscViewerASCIIPrintf(viewer,"  transpose mode: not yet set\n");CHKERRQ(ierr);
     }
+    if (svd->which == SVD_LARGEST) {
+      ierr = PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: largest\n");CHKERRQ(ierr);
+    } else {
+      ierr = PetscViewerASCIIPrintf(viewer,"  selected portion of the spectrum: smallest\n");CHKERRQ(ierr);
+    }
+    ierr = PetscViewerASCIIPrintf(viewer,"  number of singular values (nsv): %d\n",svd->nsv);CHKERRQ(ierr);
+    ierr = PetscViewerASCIIPrintf(viewer,"  number of column vectors (ncv): %d\n",svd->ncv);CHKERRQ(ierr);
+    ierr = PetscViewerASCIIPrintf(viewer,"  maximum number of iterations: %d\n",svd->max_it);
+    ierr = PetscViewerASCIIPrintf(viewer,"  tolerance: %g\n",svd->tol);CHKERRQ(ierr);
     if (svd->ops->view) {
       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
       ierr = (*svd->ops->view)(svd,viewer);CHKERRQ(ierr);

177,20 → 186,22

   svd->type_name   = PETSC_NULL;
   svd->A           = PETSC_NULL;
   svd->AT          = PETSC_NULL;
-  svd->transmode   = PETSC_DEFAULT;
+  svd->transmode   = PETSC_DECIDE;
   svd->sigma       = PETSC_NULL;
   svd->U           = PETSC_NULL;
   svd->V           = PETSC_NULL;
   svd->vec_initial = PETSC_NULL;
   svd->which       = SVD_LARGEST;
   svd->n           = 0;
-  svd->nconv       = -1;
+  svd->nconv       = 0;
   svd->nsv         = 1;
   svd->ncv         = PETSC_DECIDE;
-  svd->max_it      = PETSC_DEFAULT;
+  svd->its         = 0;
+  svd->max_it      = PETSC_DECIDE;
   svd->tol         = 1e-7;
   svd->data        = PETSC_NULL;
   svd->setupcalled = 0;
+  svd->reason      = SVD_CONVERGED_ITERATING;
   ierr = PetscPublishAll(svd);CHKERRQ(ierr);
   PetscFunctionReturn(0);

226,8 → 237,8

     ierr = (*svd->ops->destroy)(svd); CHKERRQ(ierr);
   }
-  if (svd->A) { ierr = MatDestroy(svd->A);CHKERRQ(ierr);  }
-  if (svd->AT) { ierr = MatDestroy(svd->AT);CHKERRQ(ierr);  }
+  if (svd->A) { ierr = MatDestroy(svd->A);CHKERRQ(ierr); }
+  if (svd->AT) { ierr = MatDestroy(svd->AT);CHKERRQ(ierr); }
   if (svd->n) {
     ierr = PetscFree(svd->sigma);CHKERRQ(ierr);
     for (i=0;i<svd->n;i++) {

239,7 → 250,8

     }
     ierr = PetscFree(svd->V);CHKERRQ(ierr);
   }
-  if (svd->data) { ierr = PetscFree(svd->data);CHKERRQ(ierr);  }
+  if (svd->vec_initial) { ierr = VecDestroy(svd->vec_initial);CHKERRQ(ierr); }
+  if (svd->data) { ierr = PetscFree(svd->data);CHKERRQ(ierr); }
   PetscLogObjectDestroy(svd);
   PetscHeaderDestroy(svd);

/trunk/src/svd/interface/svdopts.c

41,9 → 41,11

   PetscFunctionBegin;
   PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
   switch (mode) {
+    case PETSC_DEFAULT:
+      mode = PETSC_DECIDE;
     case SVD_TRANSPOSE_EXPLICIT:
-    case SVD_TRANSPOSE_MATMULT:
-    case PETSC_DEFAULT:
+    case SVD_TRANSPOSE_IMPLICIT:
+    case PETSC_DECIDE:
       svd->transmode = mode;
       svd->setupcalled = 0;
       break;

97,7 → 99,7

    Options Database Keys:
 +  -svd_tol <tol> - Sets the convergence tolerance
 -  -svd_max_it <maxits> - Sets the maximum number of iterations allowed
-   (use PETSC_DEFAULT to compute a value based on the operator matrix)
+   (use PETSC_DECIDE to compute an educated guess based on basis and matrix sizes)
    Notes:
    Use PETSC_IGNORE to retain the previous value of any parameter.

111,16 → 113,21

   PetscFunctionBegin;
   PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
   if (tol != PETSC_IGNORE) {
-    if (tol < 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
-    svd->tol = tol;
+    if (tol == PETSC_DEFAULT) {
+      tol = 1e-7;
+    } else {
+      if (tol < 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
+      svd->tol = tol;
+    }
   }
   if (maxits != PETSC_IGNORE) {
-    if (maxits == PETSC_DEFAULT) {
+    if (maxits == PETSC_DEFAULT || maxits == PETSC_DECIDE) {
+      svd->max_it = PETSC_DECIDE;
       svd->setupcalled = 0;
     } else {
       if (maxits < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of maxits. Must be > 0");
+      svd->max_it = maxits;
     }
-    svd->max_it = maxits;
   }
   PetscFunctionReturn(0);
 }

176,8 → 183,8

    Notes:
    Use PETSC_IGNORE to retain the previous value of any parameter.
-   Use PETSC_DEFAULT for ncv to assign a reasonably good value, which is
-   dependent on the solution method.
+   Use PETSC_DECIDE for ncv to assign a reasonably good value, which is
+   dependent on the solution method and the number of singular values required.
    Level: intermediate

191,11 → 198,14

   if (nsv != PETSC_IGNORE) {
     if (nsv<1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nsv. Must be > 0");
     svd->nsv = nsv;
-    svd->setupcalled = 0;
   }
   if (ncv != PETSC_IGNORE) {
-    if (ncv<1 && ncv != PETSC_DECIDE) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
-    svd->ncv = ncv;
+    if (ncv == PETSC_DEFAULT || ncv == PETSC_DECIDE) {
+      svd->ncv = PETSC_DECIDE;
+    } else {
+      if (ncv<1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
+      svd->ncv = ncv;
+    }
     svd->setupcalled = 0;
   }
   PetscFunctionReturn(0);

328,8 → 338,8

 {
   PetscErrorCode ierr;
   char           type[256];
-  PetscTruth     flg,flg2;
-  const char     *mode_list[2] = { "explicit", "matmult" };
+  PetscTruth     flg;
+  const char     *mode_list[2] = { "explicit", "implicit" };
   PetscInt       i,j;
   PetscReal      r;

347,24 → 357,20

   ierr = PetscOptionsName("-svd_view","Print detailed information on solver used","SVDView",0);CHKERRQ(ierr);
-  ierr = PetscOptionsEList("-svd_transpose_mode","Transpose SVD mode","SVDSetTransposeMode",mode_list,2,svd->transmode == PETSC_DEFAULT ? "default" : mode_list[svd->transmode],&i,&flg);CHKERRQ(ierr);
+  ierr = PetscOptionsEList("-svd_transpose_mode","Transpose SVD mode","SVDSetTransposeMode",mode_list,2,svd->transmode == PETSC_DECIDE ? "decide" : mode_list[svd->transmode],&i,&flg);CHKERRQ(ierr);
   if (flg) {
     ierr = SVDSetTransposeMode(svd,(SVDTransposeMode)i);CHKERRQ(ierr);
   }
   r = i = PETSC_IGNORE;
-  ierr = PetscOptionsInt("-svd_max_it","Maximum number of iterations","SVDSetTolerances",svd->max_it,&i,&flg);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-svd_tol","Tolerance","SVDSetTolerances",svd->tol,&r,&flg2);CHKERRQ(ierr);
-  if (flg || flg2) {
-    ierr = SVDSetTolerances(svd,r,i);CHKERRQ(ierr);
-  }
+  ierr = PetscOptionsInt("-svd_max_it","Maximum number of iterations","SVDSetTolerances",svd->max_it,&i,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-svd_tol","Tolerance","SVDSetTolerances",svd->tol,&r,PETSC_NULL);CHKERRQ(ierr);
+  ierr = SVDSetTolerances(svd,r,i);CHKERRQ(ierr);
   i = j = PETSC_IGNORE;
-  ierr = PetscOptionsInt("-svd_nsv","Number of singular values to compute","SVDSetDimensions",svd->ncv,&i,&flg);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-svd_ncv","Number of basis vectors","SVDSetDimensions",svd->ncv,&j,&flg2);CHKERRQ(ierr);
-  if (flg || flg2) {
-    ierr = SVDSetDimensions(svd,i,j);CHKERRQ(ierr);
-  }
+  ierr = PetscOptionsInt("-svd_nsv","Number of singular values to compute","SVDSetDimensions",svd->ncv,&i,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-svd_ncv","Number of basis vectors","SVDSetDimensions",svd->ncv,&j,PETSC_NULL);CHKERRQ(ierr);
+  ierr = SVDSetDimensions(svd,i,j);CHKERRQ(ierr);
   ierr = PetscOptionsTruthGroupBegin("-svd_largest","compute largest singular values","SVDSetWhichSingularTriplets",&flg);CHKERRQ(ierr);
   if (flg) { ierr = SVDSetWhichSingularTriplets(svd,SVD_LARGEST);CHKERRQ(ierr); }

/trunk/src/svd/interface/svdsolve.c

29,6 → 29,7

PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

if (!svd->setupcalled) { ierr = SVDSetUp(svd);CHKERRQ(ierr); }

svd->reason = SVD_CONVERGED_ITERATING;

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

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

41,6 → 42,75

 }
 #undef __FUNCT__
+#define __FUNCT__ "SVDGetIterationNumber"
+/*@
+   SVDGetIterationNumber - Gets the current iteration number. If the
+   call to SVDSolve() is complete, then it returns the number of iterations
+   carried out by the solution method.
+   Not Collective
+   Input Parameter:
+.  svd - the singular value solver context
+   Output Parameter:
+.  its - number of iterations
+   Level: intermediate
+   Notes:
+      During the i-th iteration this call returns i-1. If SVDSolve() is
+      complete, then parameter "its" contains either the iteration number at
+      which convergence was successfully reached, or failure was detected.
+      Call SVDGetConvergedReason() to determine if the solver converged or
+      failed and why.
+@*/
+PetscErrorCode SVDGetIterationNumber(SVD svd,int *its)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
+  PetscValidIntPointer(its,2);
+  *its = svd->its;
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "SVDGetConvergedReason"
+/*@C
+   SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
+   stopped.
+   Not Collective
+   Input Parameter:
+.  svd - the singular value solver context
+   Output Parameter:
+.  reason - negative value indicates diverged, positive value converged
+   (see SVDConvergedReason)
+   Possible values for reason:
++  SVD_CONVERGED_TOL - converged up to tolerance
+.  SVD_DIVERGED_ITS - required more than its to reach convergence
+-  SVD_DIVERGED_BREAKDOWN - generic breakdown in method
+   Level: intermediate
+   Notes: Can only be called after the call to SVDSolve() is complete.
+.seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason
+@*/
+PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
+  PetscValidIntPointer(reason,2);
+  *reason = svd->reason;
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
 #define __FUNCT__ "SVDGetConverged"
 /*@
    SVDGetConverged - Gets the number of converged singular values.

64,7 → 134,7

   PetscFunctionBegin;
   PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
   PetscValidIntPointer(nconv,2);
-  if (svd->nconv < 0) {
+  if (svd->reason == SVD_CONVERGED_ITERATING) {
     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSolve must be called first");
   }
   *nconv = svd->nconv;

104,7 → 174,7

   PetscFunctionBegin;
   PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
   PetscValidPointer(sigma,3);
-  if (svd->nconv < 0) {
+  if (svd->reason == SVD_CONVERGED_ITERATING) {
     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSolve must be called first");
   }
   if (i<0 || i>=svd->nconv) {

137,7 → 207,7

    Output Parameters:
 +  norm1 - the norm ||A*v-sigma*u||_2 where sigma is the
            singular value, u and v are the left and right singular vectors.
--  norm2 - the norm ||A*u-sigma*v||_2 with the same sigma, u and v
+-  norm2 - the norm ||A^T*u-sigma*v||_2 with the same sigma, u and v
    Note:
    The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).

155,7 → 225,7

   PetscFunctionBegin;
   PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
-  if (svd->nconv < 0) {
+  if (svd->reason == SVD_CONVERGED_ITERATING) {
     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSolve must be called first");
   }
   if (i<0 || i>=svd->nconv) {

163,7 → 233,7

   }
   ierr = MatGetVecs(svd->A,&v,&u);CHKERRQ(ierr);
-  ierr = SVDGetSingularTriplet(svd,i,&sigma,u,v);
+  ierr = SVDGetSingularTriplet(svd,i,&sigma,u,v);CHKERRQ(ierr);
   if (norm1) {
     ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
     ierr = MatMult(svd->A,v,x);CHKERRQ(ierr);

/trunk/src/svd/interface/svdsetup.c

85,7 → 85,6

PetscErrorCode SVDSetInitialVector(SVD svd,Vec vec)

{

PetscErrorCode ierr;

PetscReal norm;

PetscFunctionBegin;

PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

95,7 → 94,6

if (svd->vec_initial) {

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

}

ierr = VecNormalize(vec,&norm);CHKERRQ(ierr);

svd->vec_initial = vec;

PetscFunctionReturn(0);

}

155,7 → 153,6

int i;

PetscTruth flg;

PetscInt M,N;

PetscReal norm;

PetscFunctionBegin;

PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

172,20 → 169,19

     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSetOperator must be called first");
   /* determine how to build the transpose */
-  if (svd->transmode == PETSC_DEFAULT) {
+  if (svd->transmode == PETSC_DECIDE) {
     ierr = MatHasOperation(svd->A,MATOP_TRANSPOSE,&flg);CHKERRQ(ierr);
     if (flg) svd->transmode = SVD_TRANSPOSE_EXPLICIT;
-    else svd->transmode = SVD_TRANSPOSE_MATMULT;
+    else svd->transmode = SVD_TRANSPOSE_IMPLICIT;
   }
   /* build transpose matrix */
+  if (svd->AT) { ierr = MatDestroy(svd->AT);CHKERRQ(ierr); }
   switch (svd->transmode) {
     case SVD_TRANSPOSE_EXPLICIT:
-      if (svd->AT) { ierr = MatDestroy(svd->AT);CHKERRQ(ierr); }
       ierr = MatTranspose(svd->A,&svd->AT);CHKERRQ(ierr);
       break;
-    case SVD_TRANSPOSE_MATMULT:
-      if (svd->AT) { ierr = MatDestroy(svd->AT);CHKERRQ(ierr); }
+    case SVD_TRANSPOSE_IMPLICIT:
       svd->AT = PETSC_NULL;
       break;
     default:

196,7 → 192,6

if (!svd->vec_initial) {

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

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

ierr = VecNormalize(svd->vec_initial,&norm);CHKERRQ(ierr);

}

/* call specific solver setup */

/trunk/src/svd/impls/eigensolver/eigensolver.c

31,7 → 31,7

   ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr);
   eigen = (SVD_EIGENSOLVER *)svd->data;
   switch (eigen->mode) {
-    case SVDEIGENSOLVER_DIRECT:
+    case SVDEIGENSOLVER_ATA:
       ierr = MatMult(svd->A,x,eigen->x1);CHKERRQ(ierr);
       if (svd->AT) {
         ierr = MatMult(svd->AT,eigen->x1,y);CHKERRQ(ierr);

39,7 → 39,7

         ierr = MatMultTranspose(svd->A,eigen->x1,y);CHKERRQ(ierr);
       }
       break;
-    case SVDEIGENSOLVER_TRANSPOSE:
+    case SVDEIGENSOLVER_AAT:
       if (svd->AT) {
         ierr = MatMult(svd->AT,x,eigen->x1);CHKERRQ(ierr);
       } else {

106,12 → 106,12

   ierr = MatGetLocalSize(svd->A,&m,&n);CHKERRQ(ierr);
   switch (eigen->mode) {
-    case SVDEIGENSOLVER_DIRECT:
+    case SVDEIGENSOLVER_ATA:
       ierr = MatGetVecs(svd->A,PETSC_NULL,&eigen->x1);CHKERRQ(ierr);
       eigen->x2 = eigen->y1 = eigen->y2 = PETSC_NULL;
       ierr = MatCreateShell(svd->comm,n,n,PETSC_DETERMINE,PETSC_DETERMINE,svd,&eigen->mat);CHKERRQ(ierr);
       break;
-    case SVDEIGENSOLVER_TRANSPOSE:
+    case SVDEIGENSOLVER_AAT:
       ierr = MatGetVecs(svd->A,&eigen->x1,PETSC_NULL);CHKERRQ(ierr);
       eigen->x2 = eigen->y1 = eigen->y2 = PETSC_NULL;
       ierr = MatCreateShell(svd->comm,m,m,PETSC_DETERMINE,PETSC_DETERMINE,svd,&eigen->mat);CHKERRQ(ierr);

155,9 → 155,9

   ierr = EPSSetWhichEigenpairs(eigen->eps,svd->which == SVD_LARGEST ? EPS_LARGEST_REAL : EPS_SMALLEST_MAGNITUDE);CHKERRQ(ierr);
   ierr = EPSSolve(eigen->eps);CHKERRQ(ierr);
   ierr = EPSGetConverged(eigen->eps,&svd->nconv);CHKERRQ(ierr);
+  ierr = EPSGetConvergedReason(eigen->eps,(EPSConvergedReason*)&svd->reason);
   switch (eigen->mode) {
-    case SVDEIGENSOLVER_DIRECT:
+    case SVDEIGENSOLVER_ATA:
       for (i=0;i<svd->nconv;i++) {
         ierr = EPSGetEigenpair(eigen->eps,i,&sigma,PETSC_NULL,svd->V[i],PETSC_NULL);CHKERRQ(ierr);
         svd->sigma[i] = sqrt(PetscRealPart(sigma));

165,7 → 165,7

         ierr = VecScale(svd->U[i],1.0/svd->sigma[i]);CHKERRQ(ierr);
       }
       break;
-    case SVDEIGENSOLVER_TRANSPOSE:
+    case SVDEIGENSOLVER_AAT:
       for (i=0;i<svd->nconv;i++) {
         ierr = EPSGetEigenpair(eigen->eps,i,&sigma,PETSC_NULL,svd->U[i],PETSC_NULL);CHKERRQ(ierr);
         svd->sigma[i] = sqrt(PetscRealPart(sigma));

216,7 → 216,7

   PetscErrorCode  ierr;
   SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
   PetscTruth      flg;
-  const char      *mode_list[3] = { "direct" , "transpose", "cyclic" };
+  const char      *mode_list[3] = { "ata" , "aat", "cyclic" };
   PetscInt        mode;
   PetscFunctionBegin;

238,8 → 238,10

   PetscFunctionBegin;
   switch (eigen->mode) {
-    case SVDEIGENSOLVER_DIRECT:
-    case SVDEIGENSOLVER_TRANSPOSE:
+    case PETSC_DEFAULT:
+      mode = SVDEIGENSOLVER_CYCLIC;
+    case SVDEIGENSOLVER_ATA:
+    case SVDEIGENSOLVER_AAT:
     case SVDEIGENSOLVER_CYCLIC:
       eigen->mode = mode;
       svd->setupcalled = 0;

434,7 → 436,7

 {
   PetscErrorCode  ierr;
   SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
-  const char      *mode_list[3] = { "direct" , "transpose", "cyclic" };
+  const char      *mode_list[3] = { "ata" , "aat", "cyclic" };
   PetscFunctionBegin;
   ierr = PetscViewerASCIIPrintf(viewer,"eigensolver SVD mode: %s\n",mode_list[eigen->mode]);CHKERRQ(ierr);

/trunk/src/svd/impls/lanczos/gklanczos.c

8,6 → 8,7

#include "src/svd/svdimpl.h" /*I "slepcsvd.h" I*/

#include "slepcblaslapack.h"

#undef __FUNCT__

#define __FUNCT__ "SVDSetUp_LANCZOS"

19,9 → 20,9

   PetscFunctionBegin;
   ierr = MatGetSize(svd->A,&M,&N);CHKERRQ(ierr);
   if (svd->ncv == PETSC_DECIDE)
-    svd->ncv = PetscMin(PetscMin(M,N),10);
-  if (svd->max_it == PETSC_DEFAULT)
-    svd->max_it = PetscMax(PetscMax(M,N)/svd->ncv,10);
+    svd->ncv = PetscMin(PetscMin(M,N),PetscMax(2*svd->nsv,10));
+  if (svd->max_it == PETSC_DECIDE)
+    svd->max_it = PetscMax(PetscMin(M,N)/svd->ncv,100);
   PetscFunctionReturn(0);
 }

101,7 → 102,10

   ierr = VecNormalize(V[0],&norm1);CHKERRQ(ierr);
   svd->nconv = 0;
-  for (svd->its=1;svd->its<svd->max_it;svd->its++) {
+  svd->its = 0;
+  while (svd->reason == SVD_CONVERGED_ITERATING) {
+    svd->its++;
     for (i=svd->nconv;i<svd->n;i++) {
       ierr = MatMult(svd->A,V[i],U[i]);CHKERRQ(ierr);
 //      if (i>svd->nconv) { ierr = VecAXPY(U[i],-beta[i-1],U[i-1]);CHKERRQ(ierr); }

125,7 → 129,7

     ierr = PetscMemzero(Q,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);
     for (i=0;i<n;i++)
       PT[i*n+i] = Q[i*n+i] = 1.0;
-    dbdsqr_("U",&n,&n,&n,&zero,alpha+svd->nconv,beta+svd->nconv,PT,&n,Q,&n,PETSC_NULL,&n,work,&info,1);
+    LAPACKbdsqr_("U",&n,&n,&n,&zero,alpha+svd->nconv,beta+svd->nconv,PT,&n,Q,&n,PETSC_NULL,&n,work,&info,1);
     k = svd->nconv;
     for (i=svd->nconv;i<svd->n;i++) {

140,7 → 144,7

       ierr = VecMAXPY(svd->U[i],n,Q+(i-svd->nconv)*n,U+svd->nconv);CHKERRQ(ierr);
       ierr = computeres(svd,svd->sigma[i],svd->U[i],svd->V[i],&norm1,&norm2);CHKERRQ(ierr);
-      printf("[%i] sigma[%i] = %g error = %g,%g\n",svd->its,i,svd->sigma[i],norm1,norm2);
+//      printf("[%i] sigma[%i] = %g error = %g,%g\n",svd->its,i,svd->sigma[i],norm1,norm2);
       if (sqrt(norm1*norm1+norm2*norm2) < svd->tol) {
         k++;
       } else break;

150,8 → 154,11

       ierr = VecCopy(svd->U[i],U[i]);CHKERRQ(ierr);
     }
     svd->nconv = k;
-    if (svd->nconv >= svd->nsv) break;
-    ierr = VecCopy(svd->V[svd->nconv],V[svd->nconv]);CHKERRQ(ierr);
+    if (svd->its > svd->max_it) svd->reason = SVD_DIVERGED_ITS;
+    if (svd->nconv >= svd->nsv) svd->reason = SVD_CONVERGED_TOL;
+    if (svd->reason == SVD_CONVERGED_ITERATING) {
+      ierr = VecCopy(svd->V[svd->nconv],V[svd->nconv]);CHKERRQ(ierr);
+    }
   }
   ierr = VecDestroyVecs(V,n+1);CHKERRQ(ierr);

/trunk/src/svd/impls/lapack/svdlapack.c

35,7 → 35,7

   ierr = MatConvert(svd->A,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat);CHKERRQ(ierr);
   ierr = MatGetArray(mat,&pmat);CHKERRQ(ierr);
   ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
-  svd->nconv = n = PetscMin(M,N);
+  n = PetscMin(M,N);
   if (M>=N) {
      pU = PETSC_NULL;
      ierr = PetscMalloc(sizeof(PetscScalar)*N*N,&pVT);CHKERRQ(ierr);

49,9 → 49,9

   lwork = -1;
 #if defined(PETSC_USE_COMPLEX)
   ierr = PetscMalloc(sizeof(PetscReal)*(5*n*n+7*n),&rwork);CHKERRQ(ierr);
-  LAPACKgesdd("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,&qwork,&lwork,rwork,iwork,&info,1);
+  LAPACKgesdd_("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,&qwork,&lwork,rwork,iwork,&info,1);
 #else
-  LAPACKgesdd("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,&qwork,&lwork,iwork,&info,1);
+  LAPACKgesdd_("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,&qwork,&lwork,iwork,&info,1);
 #endif
   if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGESDD %d",info);
   lwork = qwork;

59,10 → 59,10

   /* computation */
 #if defined(PETSC_USE_COMPLEX)
-  LAPACKgesdd("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,work,&lwork,rwork,iwork,&info,1);
+  LAPACKgesdd_("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,work,&lwork,rwork,iwork,&info,1);
   ierr = PetscFree(rwork);CHKERRQ(ierr);
 #else
-  LAPACKgesdd("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,work,&lwork,iwork,&info,1);
+  LAPACKgesdd_("O",&M,&N,pmat,&M,svd->sigma,pU,&M,pVT,&N,work,&lwork,iwork,&info,1);
 #endif
   if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGESDD %d",info);
   ierr = PetscFree(iwork);CHKERRQ(ierr);

87,6 → 87,9

     ierr = VecRestoreArray(svd->V[i],&pv);CHKERRQ(ierr);
   }
+  svd->nconv = n;
+  svd->reason = SVD_CONVERGED_TOL;
   ierr = MatRestoreArray(mat,&pmat);CHKERRQ(ierr);
   ierr = MatDestroy(mat);CHKERRQ(ierr);
   if (M>=N) {

105,8 → 108,8

   PetscFunctionBegin;
   svd->ops->setup = SVDSetup_LAPACK;
   svd->ops->solve = SVDSolve_LAPACK;
-  if (svd->transmode == PETSC_DEFAULT)
-    svd->transmode = SVD_TRANSPOSE_MATMULT; /* don't build the transpose */
+  if (svd->transmode == PETSC_DECIDE)
+    svd->transmode = SVD_TRANSPOSE_IMPLICIT; /* don't build the transpose */
   PetscFunctionReturn(0);
 }
 EXTERN_C_END

/trunk/src/svd/svdimpl.h

36,6 → 36,7

void *data; /* placeholder for misc stuff associated

with a particular solver */

int setupcalled;

SVDConvergedReason reason;

};

EXTERN PetscErrorCode SVDRegisterAll(char *);

Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 1282 → Rev 1283