Commits

BarryFSmith committed d6f8f2f

KSPDGMRES: fix conditional variable declarations

Incorporate local BLAS/LAPACK variables with each #if !defined(missing)
case rather than trying to conditionally declare the variables at the
top of the function which is messy when several possibly missing
blas/lapack functions are used.

Reported-by: Loris Bennett <loris.bennett@fu-berlin.de>

Comments (0)

Files changed (1)

src/ksp/ksp/impls/gmres/dgmres/dgmres.c

   PetscInt       N        = dgmres->max_k+1;
   PetscInt       n        = dgmres->it+1;
   PetscReal      alpha;
-#if !defined(PETSC_MISSING_LAPACK_GETRF)
-  PetscBLASInt info;
-#endif
 
   PetscFunctionBegin;
   ierr = PetscLogEventBegin(KSP_DGMRESComputeDeflationData, ksp, 0,0,0);CHKERRQ(ierr);
 #if defined(PETSC_MISSING_LAPACK_GETRF)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&nr, &nr, TTF, &bmax, INVP, &info));
-  if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGETRF INFO=%d",(int) info);
+  {
+    PetscBLASInt info;
+    PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&nr, &nr, TTF, &bmax, INVP, &info));
+    if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGETRF INFO=%d",(int) info);
+  }
 #endif
 
   /* Save X in U and MX in MU for the next cycles and increase the size of the invariant subspace */
   PetscBLASInt   liwork;
   PetscScalar    *Ht;           /* Transpose of the Hessenberg matrix */
   PetscScalar    *t;            /* Store the result of the solution of H^T*t=h_{m+1,m}e_m */
-  PetscBLASInt   nrhs = 1;      /* Number of right hand side */
   PetscBLASInt   *ipiv;         /* Permutation vector to be used in LAPACK */
-#if !defined(PETSC_MISSING_LAPACK_HSEQR) || !defined(PETSC_MISSING_LAPACK_TRSEN)
-  PetscBLASInt ilo = 1;
-  PetscBLASInt info;
-  PetscReal    CondEig;         /* lower bound on the reciprocal condition number for the selected cluster of eigenvalues */
-  PetscReal    CondSub;         /* estimated reciprocal condition number of the specified invariant subspace. */
-  PetscBool    flag;            /* determine whether to use Ritz vectors or harmonic Ritz vectors */
-#endif
+  PetscBool      flag;            /* determine whether to use Ritz vectors or harmonic Ritz vectors */
 
   PetscFunctionBegin;
   ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr);
 #if   defined(PETSC_MISSING_LAPACK_GESV)
     SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GESV - Lapack routine is unavailable.");
 #else
-    PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&bn, &nrhs, Ht, &bn, ipiv, t, &bn, &info));
-    if (info) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB, "Error while calling the Lapack routine DGESV");
+    {
+      PetscBLASInt info;
+      PetscBLASInt nrhs = 1;
+      PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&bn, &nrhs, Ht, &bn, ipiv, t, &bn, &info));
+      if (info) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB, "Error while calling the Lapack routine DGESV");
+    }
 #endif
     /* Now form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
     for (i = 0; i < bn; i++) A[(bn-1)*bn+i] += t[i];
 #if defined(PETSC_MISSING_LAPACK_HSEQR)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"HSEQR - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S", "I", &bn, &ilo, &ihi, A, &ldA, wr, wi, Q, &ldQ, work, &lwork, &info));
-  if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XHSEQR %d",(int) info);
+  {
+    PetscBLASInt info;
+    PetscBLASInt ilo = 1;
+    PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S", "I", &bn, &ilo, &ihi, A, &ldA, wr, wi, Q, &ldQ, work, &lwork, &info));
+    if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XHSEQR %d",(int) info);
+  }
 #endif
   ierr = PetscFree(work);CHKERRQ(ierr);
 
 #if defined(PETSC_MISSING_LAPACK_TRSEN)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"TRSEN - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKtrsen",LAPACKtrsen_("B", "V", select, &bn, A, &ldA, Q, &ldQ, wr, wi, &NbrEig, &CondEig, &CondSub, work, &lwork, iwork, &liwork, &info));
-  if (info == 1) SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "UNABLE TO REORDER THE EIGENVALUES WITH THE LAPACK ROUTINE : ILL-CONDITIONED PROBLEM");
+  {
+    PetscBLASInt info;
+    PetscReal    CondEig;         /* lower bound on the reciprocal condition number for the selected cluster of eigenvalues */
+    PetscReal    CondSub;         /* estimated reciprocal condition number of the specified invariant subspace. */
+    PetscStackCallBLAS("LAPACKtrsen",LAPACKtrsen_("B", "V", select, &bn, A, &ldA, Q, &ldQ, wr, wi, &NbrEig, &CondEig, &CondSub, work, &lwork, iwork, &liwork, &info));
+    if (info == 1) SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "UNABLE TO REORDER THE EIGENVALUES WITH THE LAPACK ROUTINE : ILL-CONDITIONED PROBLEM");
+  }
 #endif
   ierr = PetscFree(select);CHKERRQ(ierr);
 
   KSP_DGMRES     *dgmres = (KSP_DGMRES*) ksp->data;
   PetscInt       i, r     = dgmres->r;
   PetscErrorCode ierr;
-  PetscBLASInt   nrhs     = 1;
   PetscReal      alpha    = 1.0;
   PetscInt       max_neig = dgmres->max_neig;
   PetscBLASInt   br,bmax;
   PetscReal      lambda = dgmres->lambdaN;
-#if !defined(PETSC_MISSING_LAPACK_GERFS)
-  PetscReal    berr, ferr;
-  PetscBLASInt info;
-#endif
 
   PetscFunctionBegin;
   ierr = PetscBLASIntCast(r,&br);CHKERRQ(ierr);
 #if defined(PETSC_MISSING_LAPACK_GETRS)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N", &br, &nrhs, TTF, &bmax, INVP, X1, &bmax, &info));
-  if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGETRS %d", (int) info);
+  {
+    PetscBLASInt info;
+    PetscBLASInt nrhs = 1;
+    PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N", &br, &nrhs, TTF, &bmax, INVP, X1, &bmax, &info));
+    if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGETRS %d", (int) info);
+  }
 #endif
   /* Iterative refinement -- is it really necessary ?? */
   if (!WORK) {
 #if defined(PETSC_MISSING_LAPACK_GERFS)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GERFS - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKgerfs",LAPACKgerfs_("N", &br, &nrhs, TT, &bmax, TTF, &bmax, INVP, X2, &bmax,X1, &bmax, &ferr, &berr, WORK, IWORK, &info));
-  if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGERFS %d", (int) info);
+  {
+    PetscBLASInt info;
+    PetscReal    berr, ferr;
+    PetscBLASInt nrhs = 1;
+    PetscStackCallBLAS("LAPACKgerfs",LAPACKgerfs_("N", &br, &nrhs, TT, &bmax, TTF, &bmax, INVP, X2, &bmax,X1, &bmax, &ferr, &berr, WORK, IWORK, &info));
+    if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGERFS %d", (int) info);
+  }
 #endif
 
   for (i = 0; i < r; i++) X2[i] =  X1[i]/lambda - X2[i];
   PetscReal    *Z;             /*  orthogonal matrix of  (right) schur vectors */
   PetscReal    *work;          /* working vector */
   PetscBLASInt lwork;          /* size of the working vector */
-  PetscBLASInt info;           /* Output info from Lapack routines */
   PetscInt     *perm;          /* Permutation vector to sort eigenvalues */
   PetscReal    *wr, *wi, *beta, *modul; /* Real and imaginary part and modul of the eigenvalues of A*/
   PetscInt     ierr;
   PetscBLASInt NbrEig = 0,nr,bm;
   PetscBLASInt *select;
   PetscBLASInt liwork, *iwork;
-#if !defined(PETSC_MISSING_LAPACK_TGSEN)
-  PetscReal    Dif[2];
-  PetscBLASInt ijob  = 2;
-  PetscBLASInt wantQ = 1, wantZ = 1;
-#endif
 
   PetscFunctionBegin;
   /* Block construction of the matrices AUU=(AU)'*U and (AU)'*AU*/
 #if defined(PETSC_MISSING_LAPACK_GGES)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GGES - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKgges",LAPACKgges_("V", "V", "N", NULL, &N, AUAU, &ldA, AUU, &ldA, &i, wr, wi, beta, Q, &N, Z, &N, work, &lwork, NULL, &info));
-  if (info) SETERRQ1 (PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGGES %d", (int) info);
+  {
+    PetscBLASInt info;
+    PetscStackCallBLAS("LAPACKgges",LAPACKgges_("V", "V", "N", NULL, &N, AUAU, &ldA, AUU, &ldA, &i, wr, wi, beta, Q, &N, Z, &N, work, &lwork, NULL, &info));
+    if (info) SETERRQ1 (PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGGES %d", (int) info);
+  }
 #endif
   for (i=0; i<N; i++) {
     if (beta[i] !=0.0) {
 #if defined(PETSC_MISSING_LAPACK_TGSEN)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"TGSEN - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKtgsen",LAPACKtgsen_(&ijob, &wantQ, &wantZ, select, &N, AUAU, &ldA, AUU, &ldA, wr, wi, beta, Q, &N, Z, &N, &NbrEig, NULL, NULL, &(Dif[0]), work, &lwork, iwork, &liwork, &info));
-  if (info == 1) SETERRQ(PetscObjectComm((PetscObject)ksp), -1, "UNABLE TO REORDER THE EIGENVALUES WITH THE LAPACK ROUTINE : ILL-CONDITIONED PROBLEM");
+  {
+    PetscBLASInt info;
+    PetscReal    Dif[2];
+    PetscBLASInt ijob  = 2;
+    PetscBLASInt wantQ = 1, wantZ = 1;
+    PetscStackCallBLAS("LAPACKtgsen",LAPACKtgsen_(&ijob, &wantQ, &wantZ, select, &N, AUAU, &ldA, AUU, &ldA, wr, wi, beta, Q, &N, Z, &N, &NbrEig, NULL, NULL, &(Dif[0]), work, &lwork, iwork, &liwork, &info));
+    if (info == 1) SETERRQ(PetscObjectComm((PetscObject)ksp), -1, "UNABLE TO REORDER THE EIGENVALUES WITH THE LAPACK ROUTINE : ILL-CONDITIONED PROBLEM");
+  }
 #endif
   ierr = PetscFree(select);CHKERRQ(ierr);
 
 #if defined(PETSC_MISSING_LAPACK_GETRF)
   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
 #else
-  PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&nr, &nr, TTF, &bm, INVP, &info));
-  if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGETRF INFO=%d",(int) info);
+  {
+    PetscBLASInt info;
+    PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&nr, &nr, TTF, &bm, INVP, &info));
+    if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGETRF INFO=%d",(int) info);
+  }
 #endif
   /* Free Memory */
   ierr = PetscFree(wr);CHKERRQ(ierr);
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.