Fixed memory leak when reading sparse matrix. Improved BiCGStab exit conditions.

BiCGStab_Iker/BiCGStab.c

@@ -180,18 +180,18 @@ void BiCGStab_compute (Compute_data *computeData, user_redist_t *user_data) {
     double DONE = 1.0, DMONE = -1.0, DZERO = 0.0;
     int n, n_dist;
     int maxiter, myId, reconfigure, rec_iter, state, flag;
-    double beta, alpha, umbral, omega, tmp;
+    double beta, alpha, umbral, omega, tmp, snrm2;
     double t3, t4;
     double reduce[2];
 
     n = size; n_dist = sizeR; maxiter = 16 * size; rec_iter = maxiter / 2; umbral = 1.0e-8;
     myId = computeData->myId;
     state = -1;
-    reconfigure = 0; rec_iter= 500;
-//	  flag = (computeData->rho == 0.0)? -1: 1;
+    reconfigure = 0; rec_iter = 500; maxiter = 1000;
+	  flag = (computeData->rho == 0.0)? -1: 1;
 
-    while ((computeData->iter < maxiter) && (computeData->tol > umbral)) {
-//    while ((computeData->iter < maxiter) && (flag == 1)) {
+//    while ((computeData->iter < maxiter) && (computeData->tol > umbral)) {
+    while ((computeData->iter < maxiter) && (flag == 1)) {
 
 #if PRECOND
         VvecDoubles (DONE, computeData->diags, computeData->p, DZERO, computeData->p_hat, n_dist);              // p_hat = D^-1 * p
@@ -211,10 +211,9 @@ void BiCGStab_compute (Compute_data *computeData, user_redist_t *user_data) {
 
         if (myId == 0) 
 #if DIRECT_ERROR
-            printf ("PD=%d %d \t %g \t %g \t %g \n", computeData->numP, computeData->iter, computeData->tol, umbral, computeData->direct_err);
+            printf ("PD=%d %d \t %20.10e \t %g \t %g \n", computeData->numP, computeData->iter, computeData->tol, umbral, computeData->direct_err);
 #else        
-        printf ("%d \t %g \n", computeData->iter, computeData->tol);
-//        printf ("%d \t %20.10e \n", computeData->iter, computeData->tol);
+        printf ("%d \t %20.10e \n", computeData->iter, computeData->tol);
 #endif // DIRECT_ERROR
         alpha = rdot (&n_dist, computeData->r0, &IONE, computeData->s, &IONE);
         MPI_Allreduce (MPI_IN_PLACE, &alpha, 1, MPI_DOUBLE, MPI_SUM, computeData->comm);
@@ -224,6 +223,15 @@ void BiCGStab_compute (Compute_data *computeData, user_redist_t *user_data) {
         rcopy (&n_dist, computeData->r, &IONE, computeData->q, &IONE);                            // q = r
         tmp = -alpha;
         raxpy (&n_dist, &tmp, computeData->s, &IONE, computeData->q, &IONE);                      // q = r - alpha * s;
+	
+        snrm2 = rnrm2 (&n_dist, computeData->s, &IONE);                               // snrm2 = norm (s)
+        MPI_Allreduce (MPI_IN_PLACE, &snrm2, 1, MPI_DOUBLE, MPI_SUM, computeData->comm);
+        if (snrm2 < umbral) {
+          raxpy (&n, &alpha, computeData->p_hat, &IONE, computeData->x, &IONE);               // x = x + alpha*p_hat;
+          if (myId == 0) printf ("snrm2 < umbral\n");
+          flag = 0;
+          break;
+	}
 
         // second spmv
 #if PRECOND
@@ -247,6 +255,11 @@ void BiCGStab_compute (Compute_data *computeData, user_redist_t *user_data) {
         MPI_Allreduce (MPI_IN_PLACE, reduce, 2, MPI_DOUBLE, MPI_SUM, computeData->comm);
 
         omega = reduce[0] / reduce[1];
+        if (omega == 0.0) {
+	  if (myId == 0) printf ("omega < 0.0\n");
+          flag = -2;
+          break;
+        }
 
         // x+1 = x + alpha * p + omega * q
         raxpy (&n_dist, &alpha, computeData->p_hat, &IONE, computeData->x, &IONE); 
@@ -264,6 +277,16 @@ void BiCGStab_compute (Compute_data *computeData, user_redist_t *user_data) {
 
         tmp = reduce[0];
         computeData->tol = sqrt (reduce[1]) / computeData->tol0;
+        if (computeData->tol < umbral) {
+	  if (myId == 0) printf ("tol < umbral\n");
+          flag = 0;
+          break;
+        }
+        if (tmp == 0.0) {
+	  if (myId == 0) printf ("rho == 0.0\n");
+          flag = -1;
+          break;
+        }
 
         // beta = (alpha / omega) * <r0, r+1> / <r0, r>
         beta = (alpha / omega) * (tmp / computeData->rho);
@@ -423,7 +446,8 @@ int main (int argc, char **argv) {
       if(mat_from_file) {
           if (myId == root) {
               // Creating the matrix
-              ReadMatrixHB (argv[1], &sym);
+//              ReadMatrixHB (argv[1], &sym);
+              CreateSparseMatrixHB (argv[1], &sym, 1);
               TransposeSparseMatrices (sym, 0, &mat, 0);
               dim = mat.dim1;
           }

BiCGStab_Iker/ScalarVectors.c

@@ -93,6 +93,54 @@ void ScaleInts (int *vint, int scal, int dim) {
     *(pi++) *= scal;
 }
 
+void GetIntFromString (char *string, int *pnum, int numC, int shft) {
+  int j = 0, num = 0, neg = 0;
+  char *pchar = string;
+
+  while ((j < numC) && ((*pchar < '0') || (*pchar > '9')) &&
+         (*pchar != '+') && (*pchar != '-')) { j++; pchar++; }
+  if (j < numC) {
+    if ((*pchar == '+') || (*pchar == '-'))
+      { neg = (*pchar == '-'); j++; pchar++; }
+    while ((j < numC) && (*pchar >= '0') && (*pchar <= '9'))
+      { num = num * 10 + (*pchar - 48); j++; pchar++; }
+  }
+  if (neg) num = -num;
+  *pnum = num + shft;
+}
+
+void GetIntsFromString (char *string, int *vec, int numN, int numC, int shft) {
+  int i, *pint = vec;
+  char *pchar = string;
+
+  for (i=0; i<numN; i++)
+    { GetIntFromString (pchar, (pint++), numC, shft); pchar += numC; }
+}
+
+void GetFormatsFromString (char *string, int *vec, int numN, int numC) {
+  int i, k = 0;
+  int *pint = vec;
+  char *pchar = string, *pch = NULL, c = ' ', c2;
+
+  for (i=0; i<numN; i++) {
+    pch = pchar;
+    while (*pch == ' ') pch++;
+    sscanf (pch, "(%i%c", pint, &c);
+    if ((c == 'P') || (c == 'p')) {
+      sscanf (pch, "(%i%c%i%c%i.%i)", &k, &c2, pint, &c, pint+1, pint+2);
+      pint += 3;
+    }
+    else if ((c == 'E') || (c == 'e') || (c == 'D') || (c == 'd') ||
+             (c == 'F') || (c == 'f') || (c == 'G') || (c == 'g')) {
+      sscanf (pch, "(%i%c%i.%i)", pint, &c, pint+1, pint+2);
+      pint += 3;
+    }
+    else
+      { sscanf (pch, "(%i%c%i)", pint, &c, pint+1); pint += 2; }
+    pchar += numC;
+  }
+}
+
 /*********************************************************************************/
 
 void CreateDoubles (double **vdbl, int dim) {
@@ -152,5 +200,59 @@ void VvecDoubles (double alfa, double *src1, double *src2, double beta, double *
     }
 }
 
+void GetDoubleFromString (char *string, double *pdbl, int numC) {
+  int j, k, exp, neg;
+  double num, frac;
+  char *pchar = string;
+
+  j = 0; exp = 0; neg = 0; num = 0.0; frac = 1.0;
+  while ((j < numC) && ((*pchar < '0') || (*pchar > '9')) &&
+         (*pchar != '+') && (*pchar != '-') && (*pchar != '.')) { j++; pchar++; }
+  if (j < numC) {
+    if ((*pchar == '+') || (*pchar == '-'))
+      { neg = (*pchar == '-'); j++; pchar++; }
+    if (j < numC) {
+      if (*pchar != '.')
+        while ((j < numC) && (*pchar >= '0') && (*pchar <= '9'))
+          { num = num * 10 + (*pchar - 48); j++; pchar++; }
+      if (j < numC) {
+        if (*pchar == '.') {
+          j++; pchar++;
+          while ((j < numC) && (*pchar >= '0') && (*pchar <= '9'))
+            { frac /= 10; num += (*pchar-48) * frac; j++; pchar++; }
+        }
+        if (neg) num = -num;
+        if (j < numC) {
+          if ((*pchar == 'e') || (*pchar == 'E') || (*pchar == 'd') || (*pchar == 'D')) {
+            neg = 0; j++; pchar++;
+            if (j < numC) {
+              if ((*pchar == '+') || (*pchar == '-'))
+                { neg = (*pchar == '-'); j++; pchar++; }
+              if (j < numC) {
+                while ((j < numC) && (*pchar >= '0') && (*pchar <= '9'))
+                  { exp = exp*10 + (*pchar-48); j++; pchar++; }
+                if (neg) exp = -exp;
+                for (k=0; k<exp; k++) num *= 10;
+                for (k=0; k>exp; k--) num /= 10;
+              }
+            }
+          }
+        }
+      } else
+        if (neg) num = -num;
+    }
+  }
+  *pdbl = num;
+}
+
+void GetDoublesFromString (char *string, double *vec, int numN, int numC) {
+  int i;
+  double *paux = vec;
+  char *pchar = string;
+
+  for (i=0; i<numN; i++)
+    { GetDoubleFromString (pchar, (paux++), numC); pchar += numC; }
+}
 
 /*********************************************************************************/
+

BiCGStab_Iker/ScalarVectors.h

@@ -30,6 +30,12 @@ extern void ComputeInvPermutation (int *perm, int *iperm, int dim);
 // Scale by scal the elements of vint, whose size is dim. 
 extern void ScaleInts (int *vint, int scal, int dim);
 
+extern void GetIntFromString (char *string, int *pnum, int numC, int shft);
+
+extern void GetIntsFromString (char *string, int *vec, int numN, int numC, int shft);
+
+extern void GetFormatsFromString (char *string, int *vec, int numN, int numC);
+
 /*********************************************************************************/
 
 extern void CreateDoubles (double **vdbl, int dim);
@@ -48,4 +54,8 @@ extern double DotDoubles (double *vdbl1, double *vdbl2, int dim);
 
 extern void VvecDoubles (double alfa, double *src1, double *src2, double beta, double *dst, int dim);
 
+extern void GetDoubleFromString (char *string, double *pdbl, int numC);
+
+extern void GetDoublesFromString (char *string, double *vec, int numN, int numC);
+
 /*********************************************************************************/

BiCGStab_Iker/SparseProduct.c

@@ -260,6 +260,79 @@ int ReadMatrixHB (char *filename, ptr_SparseMatrix p_spr) {
 	return 0;
 }
 
+#define length1 82
+#define length2 82
+
+FILE *OpenFile (char *name, char *attr) {
+  FILE *fich;
+  if ((fich = fopen (name, attr)) == NULL)
+    { printf ("File %s not exists \n", name); exit(1); }
+  return fich;
+}
+
+void ReadStringFile (FILE *file, char *string, int length) {
+  char *s = NULL;
+  if ((s = fgets (string, length, file)) == NULL)
+    { printf ("Error en lectura \n"); exit (1); }
+}
+
+void CreateSparseMatrixHB (char *nameFile, ptr_SparseMatrix spr, int FtoC) {
+  FILE *file;
+  char string[length1], *s = NULL;
+  int i, j, k = 0, shft = (FtoC)?-1:0;
+  int *vptr = NULL, *vpos = NULL;
+  double *vval = NULL;
+  int lines[5], dim[4], formats[10];
+
+  file = OpenFile (nameFile, "r");
+  ReadStringFile (file, string, length1);
+  ReadStringFile (file, string, length1);
+  GetIntsFromString (string, lines, 5, 14, 0);
+  ReadStringFile (file, string, length1);
+  GetIntsFromString ((string+14), dim, 4, 14, 0);
+
+  CreateSparseMatrix (spr, 0, dim[0], dim[1], dim[2], 0);
+  vptr = spr->vptr; vpos = spr->vpos; vval = spr->vval;
+
+  ReadStringFile (file, string, length1);
+  GetFormatsFromString (string, formats, 2, 16);
+  GetFormatsFromString ((string+32), (formats+4), 1+(lines[4] > 0), 20);
+
+  if (lines[4] > 0) ReadStringFile (file, string, length1);
+
+  j = 0;
+  for (i = 0; i < lines[1]; i++)
+    {
+      ReadStringFile (file, string, length2);
+      k = ((dim[0] + 1) - j);
+      if (k > formats[0]) k = formats[0];
+      GetIntsFromString (string, (vptr+j), k, formats[1], shft);
+      j+=formats[0];
+    }
+
+  j = 0;
+  for (i = 0; i < lines[2]; i++)
+    {
+      ReadStringFile (file, string, length2);
+      k = (dim[2] - j);
+      if (k > formats[2]) k = formats[2];
+      GetIntsFromString (string, (vpos+j), k, formats[3], shft);
+      j+=formats[2];
+    }
+
+  j = 0;
+  for (i = 0; i < lines[3]; i++)
+    {
+      ReadStringFile (file, string, length2);
+      k = (dim[2] - j);
+      if (k > formats[4]) k = formats[4];
+      GetDoublesFromString (string, (vval+j), k, formats[5]);
+      j+=formats[4];
+    }
+
+  fclose (file);
+}
+
 /*********************************************************************************/
 
 // This routine computes the product { res += spr * vec }.

BiCGStab_Iker/SparseProduct.h

@@ -44,12 +44,18 @@ extern void DesymmetrizeSparseMatrices (SparseMatrix src, int indexS, ptr_Sparse
 // This routine creates de sparse matrix dst from the matrix spr.
 // The parameters indexS and indexD indicate, respectivaly, if 0-indexing or 1-indexing is used
 // to store the sparse matrices.
-void TransposeSparseMatrices (SparseMatrix src, int indexS, ptr_SparseMatrix dst, int indexD);
+extern void TransposeSparseMatrices (SparseMatrix src, int indexS, ptr_SparseMatrix dst, int indexD);
 
 /*********************************************************************************/
 
 extern int ReadMatrixHB (char *filename, ptr_SparseMatrix p_spr);
 
+extern FILE *OpenFile (char *name, char *attr);
+
+extern void ReadStringFile (FILE *file, char *string, int length);
+
+extern void CreateSparseMatrixHB (char *nameFile, ptr_SparseMatrix spr, int FtoC);
+
 /*********************************************************************************/
 
 // This routine computes the product { res += spr * vec }.

BiCGStab_Iker/mymkl.c

 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 #include "mymkl.h"
 /**********************************************/
 
@@ -41,5 +42,15 @@ double rdot (int *n, double *x, int *incx, double *y, int *incy) {
   return aux;
 }
 
+double rnrm2(int *n, double *x, int *incx) {
+  int i, dim = *n, ix = *incx;
+  double *px = x, sum = 0.0;
+
+  for (i = 0; i < dim; i++) {
+    sum += *px * *px; px += ix;
+  }
+  return sqrt(sum);
+}
+
 /**********************************************/
 

BiCGStab_Iker/mymkl.h

@@ -10,4 +10,6 @@ void raxpy (int *n, double *alpha, double *x, int *incx, double *y, int *incy);
 
 double rdot (int *n, double *x, int *incx, double *y, int *incy);
 
+double rnrm2(int *n, double *x, int *incx);
+
 #endif