Replace max and min macros with std functions when compiling with c++.

Add new make action to instal only the static library

makefile

@@ -71,6 +71,15 @@ python_install: python
 R_install:
 	@$(MAKE) -C R install
 
+install_static: lib
+	install -d $(includedir)
+	cd include && install -m 644 primme_eigs_f77.h primme_eigs_f90.inc primme_eigs.h  \
+	        primme_f77.h primme_f90.inc primme.h primme_svds_f77.h  \
+	        primme_svds_f90.inc primme_svds.h \
+		$(includedir)
+	install -d $(libdir)
+	install -m 644 lib/$(LIBRARY) $(libdir)
+
 install: solib
 	install -d $(includedir)
 	cd include && install -m 644 primme_eigs_f77.h primme_eigs_f90.inc primme_eigs.h  \

src/eigs/auxiliary_eigs.c

@@ -511,7 +511,7 @@ int machineEpsMatrix_Sprimme(double *eps, primme_context ctx) {
             primme->massMatrixMatvec_type, &eps_massmatvec));
    }
 
-   *eps = max(max(MACHINE_EPSILON, eps_matvec), eps_massmatvec);
+   *eps = max(max((double)MACHINE_EPSILON, eps_matvec), eps_massmatvec);
 
    return 0;
 }
@@ -546,7 +546,7 @@ int machineEpsOrth_Sprimme(double *eps, primme_context ctx) {
       }
    }
 
-   *eps = max(max(MACHINE_EPSILON, eps_globalsum), eps_broadcast);
+   *eps = max(max((double)MACHINE_EPSILON, eps_globalsum), eps_broadcast);
 
    return 0;
 }

src/eigs/auxiliary_eigs_normal.c

@@ -74,8 +74,8 @@ int Num_compute_residuals_Sprimme(PRIMME_INT m, int n, HEVAL *eval,
 #ifdef USE_HOST
    int j;
    for (j = 0; j < n; j++) {
-      int k, M = min(m, PRIMME_BLOCK_SIZE);
-      for (k = 0; k < m; k += M, M = min(M, m - k)) {
+      int k, M = min(m, (PRIMME_INT)PRIMME_BLOCK_SIZE);
+      for (k = 0; k < m; k += M, M = min((PRIMME_INT)M, m - k)) {
          CHKERR(Num_copy_Sprimme(
                M, &Ax[ldAx * j + k], 1, &r[ldr * j + k], 1, ctx));
          CHKERR(Num_axpy_Sprimme(
@@ -172,7 +172,8 @@ int Num_update_VWXR_Sprimme(SCALAR *V, SCALAR *W, SCALAR *BV, PRIMME_INT mV,
 
    PRIMME_INT i;     /* Loop variables */
    int j;            /* Loop variables */
-   int m=min(PRIMME_BLOCK_SIZE, mV);   /* Number of rows in the cache */
+   int m = min(
+         (PRIMME_INT)PRIMME_BLOCK_SIZE, mV); /* Number of rows in the cache */
    int nXb, nXe, nYb, nYe, nBXb, nBXe, ldG0=0, ldH0=0;
    PRIMME_INT ldX, ldY, ldBX;
    SCALAR *X, *Y, *BX;
@@ -250,7 +251,7 @@ int Num_update_VWXR_Sprimme(SCALAR *V, SCALAR *W, SCALAR *BV, PRIMME_INT mV,
    if (H) CHKERR(Num_zero_matrix_SHprimme(H0, nH, nH, ldH0, ctx));
    if (xnorms) for (i=nxb; i<nxe; i++) xnorms[i-nxb] = 0.0;
 
-   for (i=0; i < mV; i+=m, m=min(m,mV-i)) {
+   for (i = 0; i < mV; i += m, m = min((PRIMME_INT)m, mV - i)) {
       /* X = V*h(nXb:nXe-1) */
       CHKERR(Num_gemm_dhd_Sprimme("N", "N", m, nXe-nXb, nV, 1.0,
          &V[i], ldV, &h[nXb*ldh], ldh, 0.0, X, ldX, ctx));
@@ -361,9 +362,9 @@ int Num_update_VWXR_Sprimme(SCALAR *V, SCALAR *W, SCALAR *BV, PRIMME_INT mV,
       if (xnorms) for (i=nxb; i<nxe; i++) tmp[j++] = xnorms[i-nxb];
       if (j) CHKERR(globalSum_RHprimme(tmp, j, ctx));
       j = 0;
-      if (Rnorms) for (i=nRb; i<nRe; i++) Rnorms[i-nRb] = sqrt(tmp[j++]);
-      if (rnorms) for (i=nrb; i<nre; i++) rnorms[i-nrb] = sqrt(tmp[j++]);
-      if (xnorms) for (i=nxb; i<nxe; i++) xnorms[i-nxb] = sqrt(tmp[j++]);
+      if (Rnorms) for (i=nRb; i<nRe; i++, j++) Rnorms[i-nRb] = sqrt(tmp[j]);
+      if (rnorms) for (i=nrb; i<nre; i++, j++) rnorms[i-nrb] = sqrt(tmp[j]);
+      if (xnorms) for (i=nxb; i<nxe; i++, j++) xnorms[i-nxb] = sqrt(tmp[j]);
       CHKERR(Num_free_RHprimme(tmp, ctx));
    }
    else {

src/eigs/correction.c

@@ -279,13 +279,15 @@ int solve_correction_Sprimme(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
 
             if (primme->target == primme_smallest) {
                blockOfShifts[blockIndex] = eval - robustShift;
-               if (sortedIndex > 0) blockOfShifts[blockIndex] = 
-                  max(blockOfShifts[blockIndex], sortedRitzVals[sortedIndex-1]);
+               if (sortedIndex > 0)
+                  blockOfShifts[blockIndex] = max(blockOfShifts[blockIndex],
+                        (double)sortedRitzVals[sortedIndex - 1]);
             }
             else {
                blockOfShifts[blockIndex] = eval + robustShift;
-               if (sortedIndex > 0) blockOfShifts[blockIndex] = 
-                  min(blockOfShifts[blockIndex], sortedRitzVals[sortedIndex-1]);
+               if (sortedIndex > 0)
+                  blockOfShifts[blockIndex] = min(blockOfShifts[blockIndex],
+                        (double)sortedRitzVals[sortedIndex - 1]);
             } /* robust shifting */
 
          }  /* for loop */
@@ -584,12 +586,13 @@ STATIC HREAL computeRobustShift(int blockIndex, double resNorm,
    /* in The Symmetric Eigenvalue Problem by B.N. Parlett.                 */
 
    if (gap > resNorm) {
-      epsilon = min(
-            delta, min(resNorm * resNorm * primme->stats.estimateInvBNorm / gap,
-                         lowerGap));
+      epsilon = min((double)delta,
+            min(resNorm * resNorm * primme->stats.estimateInvBNorm / gap,
+                  (double)lowerGap));
    }
    else {
-      epsilon = min(resNorm * sqrt(primme->stats.estimateInvBNorm), lowerGap);
+      epsilon = min(
+            resNorm * sqrt(primme->stats.estimateInvBNorm), (double)lowerGap);
    }
 
    *approxOlsenShift = min(delta, epsilon);

src/eigs/init.c

@@ -177,7 +177,8 @@ int init_basis_Sprimme(SCALAR *V, PRIMME_INT nLocal, PRIMME_INT ldV, SCALAR *W,
    } else {
       initSize = min(primme->minRestartSize, primme->initSize);
    }
-   initSize = max(0, min(primme->n - primme->numOrthoConst, initSize));
+   initSize = max((PRIMME_INT)0,
+         min(primme->n - primme->numOrthoConst, (PRIMME_INT)initSize));
    *numGuesses = primme->initSize - initSize;
    *nextGuess = primme->numOrthoConst + initSize;
 
@@ -195,7 +196,8 @@ int init_basis_Sprimme(SCALAR *V, PRIMME_INT nLocal, PRIMME_INT ldV, SCALAR *W,
       break;
    default: assert(0);
    }
-   random = max(0, min(primme->n - primme->numOrthoConst - initSize, random));
+   random = max((PRIMME_INT)0,
+         min(primme->n - primme->numOrthoConst - initSize, (PRIMME_INT)random));
    for (i = 0; i < random; i++) {
       Num_larnv_Sprimme(
             2, primme->iseed, nLocal, &V[ldV * (initSize + i)], ctx);
@@ -210,8 +212,8 @@ int init_basis_Sprimme(SCALAR *V, PRIMME_INT nLocal, PRIMME_INT ldV, SCALAR *W,
    CHKERR(matrixMatvec_Sprimme(V, nLocal, ldV, W, ldW, 0, *basisSize, ctx));
 
    if (primme->initBasisMode == primme_init_krylov) {
-      int minRestartSize =
-            min(primme->minRestartSize, primme->n - primme->numOrthoConst);
+      int minRestartSize = min((PRIMME_INT)primme->minRestartSize,
+            primme->n - primme->numOrthoConst);
       CHKERR(init_block_krylov(V, nLocal, ldV, W, ldW, BV, ldBV, *basisSize,
             minRestartSize - 1, evecs, ldevecs, primme->numOrthoConst, VtBV,
             ldVtBV, fVtBV, ldfVtBV, maxRank, ctx));

src/eigs/main_iter.c

@@ -1015,9 +1015,9 @@ int main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
 
             /* Limit basisSize to the matrix dimension */
 
-            availableBlockSize =
-                  max(0, min(availableBlockSize,
-                               primme->n - numLocked - primme->numOrthoConst));
+            availableBlockSize = max((PRIMME_INT)0,
+                  min((PRIMME_INT)availableBlockSize,
+                        primme->n - numLocked - primme->numOrthoConst));
 
             // Calling prepare_candidates before restarting is compulsory when:
             // - the order of the Ritz pairs depends on the residual vector
@@ -1137,10 +1137,11 @@ int main_iter_Sprimme(HEVAL *evals, SCALAR *evecs, PRIMME_INT ldevecs,
 
             /* Don't increase basis size beyond matrix dimension */
 
-            numNew = max(0,
-                  min(basisSize+numNew+primme->numOrthoConst+numLocked,
-                     primme->n)
-                  - primme->numOrthoConst - numLocked - basisSize);
+            numNew = max((PRIMME_INT)0,
+                  min((PRIMME_INT)basisSize + numNew + primme->numOrthoConst +
+                              numLocked,
+                        primme->n) -
+                        primme->numOrthoConst - numLocked - basisSize);
 
             CHKERR(Num_copy_matrix_Sprimme(&evecs[nextGuess * ldevecs],
                   primme->nLocal, numNew, ldevecs, &V[basisSize * ldV], ldV,
@@ -1550,7 +1551,7 @@ STATIC int prepare_candidates(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
 
    if (blockNormsSize > 0) {
       for (*smallestResNorm = HUGE_VAL, i=0; i < blockNormsSize; i++) {
-         *smallestResNorm = min(*smallestResNorm, blockNorms[i]);
+         *smallestResNorm = min(*smallestResNorm, (double)blockNorms[i]);
       }
    }
 
@@ -1582,7 +1583,7 @@ STATIC int prepare_candidates(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
             evals0[i] = evals[map[i]];
             resNorms0[i] = resNorms[map[i]];
             primme->stats.maxConvTol =
-                  max(primme->stats.maxConvTol, resNorms0[i]);
+                  max(primme->stats.maxConvTol, (double)resNorms0[i]);
          } else {
             evals0[i] = HUGE_VAL;
             resNorms0[i] = HUGE_VAL;
@@ -1642,7 +1643,7 @@ STATIC int prepare_candidates(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
                resNorms[iev[blki]] = blockNorms[blki];
                if (flagsBlock[i] == CONVERGED)
                   primme->stats.maxConvTol =
-                        max(primme->stats.maxConvTol, blockNorms[blki]);
+                        max(primme->stats.maxConvTol, (double)blockNorms[blki]);
             }
 
             /* Count the new solution */
@@ -1673,7 +1674,7 @@ STATIC int prepare_candidates(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
             if (*blockSize == 0) {
                *smallestResNorm = HUGE_VAL;
             }
-            *smallestResNorm = min(*smallestResNorm, blockNorms[blki]);
+            *smallestResNorm = min(*smallestResNorm, (double)blockNorms[blki]);
 
             blockNorms[*blockSize] = blockNorms[blki];
             iev[*blockSize] = iev[blki];
@@ -1762,7 +1763,8 @@ STATIC int prepare_candidates(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
       /* Don't trust residual norm smaller than the error in the residual norm */
 
       for (i = *blockSize; i < blockNormsSize; i++) {
-         blockNorms[i] = max(blockNorms[i], primme->stats.estimateResidualError);
+         blockNorms[i] =
+               max((double)blockNorms[i], primme->stats.estimateResidualError);
       }
 
 
@@ -1772,8 +1774,8 @@ STATIC int prepare_candidates(SCALAR *V, PRIMME_INT ldV, SCALAR *W,
          for (i=0; i<blockNormsSize; i++) {
             primme->stats.estimateBNorm = max(
                   primme->stats.estimateBNorm, 1.0 / XNorms[*blockSize + i]);
-            primme->stats.estimateInvBNorm = max(
-                  primme->stats.estimateInvBNorm, XNorms[*blockSize + i]);
+            primme->stats.estimateInvBNorm = max(primme->stats.estimateInvBNorm,
+                  (double)XNorms[*blockSize + i]);
          }
       }
    }

src/eigs/ortho.c

@@ -694,7 +694,8 @@ STATIC int Bortho_block_gen_Sprimme(SCALAR *V, PRIMME_INT ldV, HSCALAR *VLtBVL,
          /* N(i) = ||Xc(:,i)||; normXc = max ||N(i)|| */
 
          for (i=0; i<b2-b1; i++) {
-            N[i] = sqrt(max(ABS(C[(b2-b1)*i+i]), eps_orth));
+            N[i] = sqrt(
+                  max((double)ABS(C[(b2 - b1) * i + i]), (double)eps_orth));
          }
          for (i = 0; i < b2 - b1; i++)
             for (j = 0; j <= i; j++) C[(b2 - b1) * i + j] /= N[i] * N[j];
@@ -706,7 +707,7 @@ STATIC int Bortho_block_gen_Sprimme(SCALAR *V, PRIMME_INT ldV, HSCALAR *VLtBVL,
          /* D = sqrt(D) */
 
          for (i=0; i<b2-b1; i++) {
-            D[i] = sqrt(max(D[i], eps_orth * (b2 - b1)));
+            D[i] = sqrt(max((double)D[i], eps_orth * (b2 - b1)));
          }
       }
 
@@ -854,11 +855,12 @@ int ortho_single_iteration_Sprimme(SCALAR *Q, int nQ, PRIMME_INT ldQ,
             "C", "N", nQ, nX, mQ, 1.0, Q, ldQ, X, ldX, 0.0, y, nQ, ctx));
    }
    else {
-      int m=min(M, mQ);
+      int m = min((PRIMME_INT)M, mQ);
       SCALAR *X0;
       CHKERR(Num_malloc_Sprimme(m*nX, &X0, ctx));
       CHKERR(Num_zero_matrix_SHprimme(y, nQ, nX, nQ, ctx));
-      for (i=0, m=min(M,mQ); i < mQ; i+=m, m=min(m,mQ-i)) {
+      for (i = 0, m = min((PRIMME_INT)M, mQ); i < mQ;
+            i += m, m = min((PRIMME_INT)m, mQ - i)) {
          CHKERR(Num_copy_matrix_columns_Sprimme(
                &X[i], m, inX, nX, ldX, X0, NULL, m, ctx));
          CHKERR(Num_gemm_ddh_Sprimme(
@@ -894,12 +896,12 @@ int ortho_single_iteration_Sprimme(SCALAR *Q, int nQ, PRIMME_INT ldQ,
    /* X = X - BQ*(QtBQ\y); norms(i) = norm(X(i)) */
 
    SCALAR *X0 = NULL;
-   int m=min(M, mQ);
+   int m = min((PRIMME_INT)M, mQ);
    if (inX) {
       CHKERR(Num_malloc_Sprimme(m*nX, &X0, ctx));
    }
    if (norms) for (i=0; i<nX; i++) norms[i] = 0.0;
-   for (i=0; i < mQ; i+=m, m=min(m,mQ-i)) {
+   for (i = 0; i < mQ; i += m, m = min((PRIMME_INT)m, mQ - i)) {
       if (inX) {
          CHKERR(Num_copy_matrix_columns_Sprimme(
                &X[i], m, inX, nX, ldX, X0, NULL, m, ctx));
@@ -975,13 +977,15 @@ STATIC int Num_ortho_kernel(SCALAR *Q, PRIMME_INT M, int nQ, PRIMME_INT ldQ,
 
    if (D && Y) {
       if (nX == 1) {
-         m = min(PRIMME_BLOCK_SIZE, M);   /* Number of rows in the cache */
+         m = min((PRIMME_INT)PRIMME_BLOCK_SIZE,
+               M); /* Number of rows in the cache */
          if (!Yortho) {
             Y[0] = (HSCALAR)1.0/Y[0];
             Yortho = 1;
          }
       } else {
-         m = min(PRIMME_BLOCK_SIZE, M);   /* Number of rows in the cache */
+         m = min((PRIMME_INT)PRIMME_BLOCK_SIZE,
+               M); /* Number of rows in the cache */
       }
       if (Yortho) {
          CHKERR(Num_malloc_Sprimme(m * nX, &Xo, ctx));
@@ -1167,7 +1171,7 @@ STATIC int decomposition(HSCALAR *H, int n, int ldH, HSCALAR *Y, int ldY,
 STATIC int rank_estimation(HSCALAR *V, int n0, int n1, int n, int ldV) {
 
    int i, j;
-   HREAL threashold = max(.8 / n, MACHINE_EPSILON);
+   HREAL threashold = max((HREAL).8 / n, MACHINE_EPSILON);
    for(i=n0; i<n1; i++) {
       HREAL Vii = REAL_PART(V[i * ldV + i]);
       if (!ISFINITE(Vii) || Vii <= (HREAL)0.0) break;
@@ -1215,7 +1219,7 @@ int orthogonality_error_Sprimme(
                              CONJ(VtV[i * ldVtV + i] - (HSCALAR)1.0));
          for (j = i + 1; j < n; j++)
             norm_i += REAL_PART(VtV[j * ldVtV + i] * CONJ(VtV[j * ldVtV + i]));
-         norm = max(norm, sqrt(norm_i));
+         norm = max(norm, (HREAL)sqrt(norm_i));
       }
    }
    CHKERR(broadcast_RHprimme(&norm, 1, ctx));

src/eigs/primme_c.c

@@ -565,11 +565,11 @@ STATIC void convTestFunAbsolute(double *eval, void *evec, double *rNorm,
    (void)evec; /* unused parameter */
 
    if (primme->massMatrixMatvec == NULL) {
-      *isConv = *rNorm < max(primme->eps, MACHINE_EPSILON * 2) *
+      *isConv = *rNorm < max(primme->eps, (double)MACHINE_EPSILON * 2) *
                                problemNorm_Sprimme(0, primme);
    }
    else {
-      *isConv = *rNorm < max(primme->eps, MACHINE_EPSILON) *
+      *isConv = *rNorm < max(primme->eps, (double)MACHINE_EPSILON) *
                                problemNorm_Sprimme(0, primme);
    }
    *ierr = 0;

src/eigs/primme_interface.c

@@ -563,14 +563,18 @@ void primme_set_defaults(primme_params *primme) {
    if (primme->maxBasisSize == 0) {
       if (primme->target==primme_smallest || primme->target==primme_largest)
          primme->maxBasisSize = min(primme->n - primme->numOrthoConst,
-               max(max(15, 4 * primme->maxBlockSize +
-                                 primme->restartingParams.maxPrevRetain),
+               (PRIMME_INT)max(
+                     (double)max(
+                           15, 4 * primme->maxBlockSize +
+                                     primme->restartingParams.maxPrevRetain),
                      2.5 * primme->minRestartSize +
                            primme->restartingParams.maxPrevRetain));
       else
          primme->maxBasisSize = min(primme->n - primme->numOrthoConst,
-               max(max(35, 5 * primme->maxBlockSize +
-                                 primme->restartingParams.maxPrevRetain),
+               (PRIMME_INT)max(
+                     (double)max(
+                           35, 5 * primme->maxBlockSize +
+                                     primme->restartingParams.maxPrevRetain),
                      1.7 * primme->minRestartSize +
                            primme->restartingParams.maxPrevRetain));
    }
@@ -588,22 +592,22 @@ void primme_set_defaults(primme_params *primme) {
       /* restart=basis-block*ceil((basis-restart-prevRetain)/block)-prevRetain*/
       if (primme->maxBlockSize > 1) {
          if (primme->restartingParams.maxPrevRetain > 0)
-            primme->minRestartSize =
-                  max(1, primme->maxBasisSize -
-                               primme->maxBlockSize *
-                                     (1 + (primme->maxBasisSize -
-                                                primme->minRestartSize - 1 -
-                                                primme->restartingParams
-                                                      .maxPrevRetain) /
-                                                 (double)primme->maxBlockSize) -
-                               primme->restartingParams.maxPrevRetain);
+            primme->minRestartSize = max(
+                  1., primme->maxBasisSize -
+                            primme->maxBlockSize *
+                                  (1 + (primme->maxBasisSize -
+                                             primme->minRestartSize - 1 -
+                                             primme->restartingParams
+                                                   .maxPrevRetain) /
+                                              (double)primme->maxBlockSize) -
+                            primme->restartingParams.maxPrevRetain);
          else
-            primme->minRestartSize =
-                  max(1, primme->maxBasisSize -
-                               primme->maxBlockSize *
-                                     (1 + (primme->maxBasisSize -
-                                                primme->minRestartSize - 1) /
-                                                 (double)primme->maxBlockSize));
+            primme->minRestartSize = max(
+                  1., primme->maxBasisSize -
+                            primme->maxBlockSize *
+                                  (1 + (primme->maxBasisSize -
+                                             primme->minRestartSize - 1) /
+                                              (double)primme->maxBlockSize));
       }
    }