JWS-112 Bumping version of ClustalO (src, binaries and windows) to version 1.2.4.

[jabaws.git] / binaries / src / clustalo / src / clustal / symmatrix.c

diff --git a/binaries/src/clustalo/src/clustal/symmatrix.c b/binaries/src/clustalo/src/clustal/symmatrix.c
index 169ee91..b37d883 100644 (file)
--- a/binaries/src/clustalo/src/clustal/symmatrix.c
+++ b/binaries/src/clustalo/src/clustal/symmatrix.c
@@ -15,7 +15,7 @@
  ********************************************************************/
 
 /*
- *  RCS $Id: symmatrix.c 230 2011-04-09 15:37:50Z andreas $
+ *  RCS $Id: symmatrix.c 303 2016-06-13 13:37:47Z fabian $
  *
  *
  * Functions for symmetric (square) matrices including diagonal.
@@ -35,8 +35,10 @@
 #include <ctype.h>
 #include <string.h>
 #include <assert.h>
+#include <stdbool.h>
 #include "symmatrix.h"
 
+#include "util.h"
 
 #if 0
 #define DEBUG
@@ -75,6 +77,10 @@ NewSymMatrix(symmatrix_t **symmat, int nrows, int ncols)
         fprintf(stderr, "Couldn't allocate memory (%s|%s)\n",
                 __FILE__, __FUNCTION__);
         return -1;
+    } 
+    else {
+        (*symmat)->nrows = nrows;
+        (*symmat)->ncols = ncols;
     }
     
     (*symmat)->data = (double **) malloc(nrows * sizeof(double *));
@@ -85,6 +91,38 @@ NewSymMatrix(symmatrix_t **symmat, int nrows, int ncols)
         *symmat = NULL;
         return -1;
     }
+
+#ifdef HAVE_OPENMP
+    /* one malloc for the full matrix data structure
+       assumes ncols >= nrows (asserted above) */
+    double *mat_data;
+    int elements_to_date=0;
+    if ((mat_data = (double *)malloc(((nrows+1)*nrows/2 + (ncols-nrows)*nrows) * sizeof(double))) == NULL) {
+        fprintf(stderr, "Couldn't allocate MPI memory (%s|%s)\n", __FILE__, __FUNCTION__); 
+        free((*symmat)->data);
+        free(*symmat);
+        *symmat = NULL;
+        return -1;
+    }
+    for (i=0; i<nrows; i++) {
+        (*symmat)->data[i] = &mat_data[elements_to_date];
+        elements_to_date += ncols-i;
+#ifdef TRACE
+        fprintf(stderr, "DEBUG(%s|%s():%d) initialising symmat with the number of the beast\n",
+                __FILE__, __FUNCTION__, __LINE__);
+#endif
+        {
+            int j;
+            for (j=0; j<ncols-i; j++) {
+#ifdef TRACE
+                (*symmat)->data[i][j] = -666.0;
+#else
+                (*symmat)->data[i][j] = 0.0;
+#endif
+            }
+        }
+    }
+#else  /* not HAVE_OPENMP */
     for (i=0; i<nrows; i++) {
         (*symmat)->data[i] = (double *) calloc((ncols-i), sizeof(double));
         if (NULL == (*symmat)->data[i]) {
@@ -110,10 +148,8 @@ NewSymMatrix(symmatrix_t **symmat, int nrows, int ncols)
         }
 #endif
     }
+#endif /* HAVE_OPENMP */
     
-    (*symmat)->nrows = nrows;
-    (*symmat)->ncols = ncols;
-
     return 0;
 }
 /***   end: new_symmatrix   ***/
@@ -240,12 +276,18 @@ SymMatrixGetValueP(double **val, symmatrix_t *symmat,
 void
 FreeSymMatrix(symmatrix_t **symmat)
 {
+#ifndef HAVE_OPENMP
     int i;
+#endif
     if (NULL != (*symmat)) {
         if (NULL != (*symmat)->data) {
+#ifdef HAVE_OPENMP
+            free((*symmat)->data[0]);
+#else
             for (i=0; i<(*symmat)->nrows; i++) {
                 free((*symmat)->data[i]);
             }
+#endif
             free((*symmat)->data);
         }
     }
@@ -271,7 +313,7 @@ FreeSymMatrix(symmatrix_t **symmat)
  *
  */
 void
-SymMatrixPrint(symmatrix_t *symmat, char **labels,  const char *path)
+SymMatrixPrint(symmatrix_t *symmat, char **labels,  const char *path, bool bPercID)
 {
     FILE *fp = NULL;
     int i, j; /* aux */
@@ -315,8 +357,15 @@ SymMatrixPrint(symmatrix_t *symmat, char **labels,  const char *path)
         /* we are lazy and do it like fastdist: write all distances
          *  for one seq to one line
          */
-        for (j=0; j<symmat->ncols; j++) {
-            fprintf(fp, " %f", SymMatrixGetValue(symmat, i, j));
+        if (TRUE == bPercID){
+            for (j=0; j<symmat->ncols; j++) {
+                fprintf(fp, " %f", 100*(1.00-SymMatrixGetValue(symmat, i, j)));
+            }
+        }
+        else{
+            for (j=0; j<symmat->ncols; j++) {
+                fprintf(fp, " %f", SymMatrixGetValue(symmat, i, j));
+            }
         }
         fprintf(fp, "%s", "\n");
     }
@@ -336,6 +385,8 @@ SymMatrixPrint(symmatrix_t *symmat, char **labels,  const char *path)
  *
  * @param[in] pcFileIn
  * distance matrix filename 
+ * @param[in] prMSeq
+ * optional mseq pointer. only used for consistency checking of names/labels 
  * @param[out] prSymMat_p
  * the symmatrix_t. will be allocated here.
  * @return:
@@ -345,7 +396,7 @@ SymMatrixPrint(symmatrix_t *symmat, char **labels,  const char *path)
  * FIXME untested code
  */
 int
-SymMatrixRead(char *pcFileIn, symmatrix_t **prSymMat_p)
+SymMatrixRead(char *pcFileIn, symmatrix_t **prSymMat_p, mseq_t *prMSeq)
 {    
     FILE *prFilePointer;
     char *buf;
@@ -359,10 +410,7 @@ SymMatrixRead(char *pcFileIn, symmatrix_t **prSymMat_p)
 
     assert(NULL!=pcFileIn);
 
-    /* FIXME: test correctness and implement label checking */
-    fprintf(stderr, "WARNING: Reading of distance matrix from file not thoroughly tested!\n");
-    fprintf(stderr, "WARNING: Assuming same order of sequences in sequence file and distance matrix file (matching of labels not implemented)\n");
-
+    /* could pass down prMSeq for label checking if non NULL */
     if (NULL == (buf = (char *) malloc(MAX_BUF_SIZE * sizeof(char)))) {
         fprintf(stderr, "ERROR: couldn't allocate memory at %s:%s:%d\n",
                 __FILE__, __FUNCTION__, __LINE__);
@@ -445,7 +493,12 @@ SymMatrixRead(char *pcFileIn, symmatrix_t **prSymMat_p)
             while (isspace(szToken[strlen(szToken)-1])) {
                 szToken[strlen(szToken)-1]='\0';
             }
-            /* FIXME save label? */
+            if (strcmp(szToken, prMSeq->sqinfo[iNParsedSeq-1].name) != 0) {
+                fprintf(stderr, "Sequence ordering in mseq and distmat differ (expected %s and got %s from distmat %s)n", 
+                        prMSeq->sqinfo[iNParsedSeq-1].name, szToken, pcFileIn);
+                iRetCode = -1;
+                goto closefile_and_freebuf;                
+            }
             szToken = strtok(NULL, " \t");
         }
         /* from here on it's only parsing of distances */