text: /****************************************************************************
   ****************************************************************************
   *                                                                          *
   *  Copyright (C) 2018  Genome Research Ltd.                                *
   *                                                                          *
   *  Author: Zemin Ning (zn1@sanger.ac.uk)                                   *
   *                                                                          *
   *  This file is part of ScaffHiC pipeline.                                 *
   *                                                                          *
   *  ScaffHiC is a free software: you can redistribute it and/or modify it   *
   *  under the terms of the GNU General Public License as published by the   *
   *  Free Software Foundation, either version 3 of the License, or (at your  *
   *  option) any later version.                                              *
   *                                                                          *
   *  This program is distributed in the hope that it will be useful, but     *
   *  WITHOUT ANY WARRANTY; without even the implied warranty of              *
   *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU        *
   *  General Public License for more details.                                *
   *                                                                          *
   *  You should have received a copy of the GNU General Public License along *
   *  with this program.  If not, see <http://www.gnu.org/licenses/>.         *
   *                                                                          *
   ****************************************************************************
   ****************************************************************************/
  /****************************************************************************/
  
  
  #include <math.h>
  #include <values.h>
  #include <stdio.h>
  #include <netinet/in.h>
  #include <stdlib.h>
  #include <dirent.h>
  #include <string.h>
  #include <ctype.h>
  #include <unistd.h>
  #include <sys/types.h> 
  #include <sys/socket.h>
  #include <sys/wait.h>
  #include <sys/signal.h>
  #include <errno.h>
  #include "fasta.h"
  
  #define MAXLINE 4096
  #define ENDS_EXTRA 0
  #define PADCHAR '-'
  #define MAX_N_BRG 50000 
  #define MAX_N_ROW 40000 
  #define Max_N_NameBase 60
  #define Max_N_Pair 100
  
  /* SSAS default parameters   */
  static int IMOD=0;
  static int num_reads=0;
  static int n_trim = 100;
  static int cons_flag = 0;
  static int min_length = 20;
  static int N_outread = 40000;
  static char strain_name[100];
  
  typedef struct
  {
         int foffset;
         int fsindex;
  } SIO;
  
  fasta *sub;
  
  static char rc_char[500];
  
  int Reverse_Complement_Contig(char c_array[],int num_len)
  {
          int i,len;
          char *tp,*dp;
  
          len=num_len;
          dp=rc_char;
          tp = c_array+len;
          for(i=len;--i>=0;)
          {
                  int tmp = *--tp;
                  if     (tmp == 't') *dp++ = 'a';
                  else if(tmp == 'g') *dp++ = 'c';
                  else if(tmp == 'c') *dp++ = 'g';
                  else if(tmp == 'a') *dp++ = 't';
                  else                *dp++ = tmp;
          }
          return(0);
  }
  
  int main(int argc, char **argv)
  {
      int i,nSeq,args;
      char **cmatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch);
      fasta *seq; 
      void ArraySort_Mix(int n, B64_long *arr, int *brr);
      void SNP_Consensus(char **argv,int args,int nSeq);
      void Read_Pairs(char **argv,int args,int nLib,int nSeq);
      void Align_Process(char **argv,int args,int argc,int nRead);
  
      seq=NULL;
  //    fflush(stdout);
  //    system("ps aux | grep eulerSNP; date");
      memset(strain_name,'\0',100);
      if(argc < 2)
      {
        printf("USAGE\n");
        printf("\n");
        printf(" %s <-split 4000> <output fastq head> <input fastq files>\n",argv[0]);
        printf("-split 4000 means cut the file into small patches each with 4000 sequences\n");
        printf("\n");
        exit(1);
      }
  
      nSeq=0;
      args=1;
      for(i=1;i<argc;i++)
      {
         if(!strcmp(argv[i],"-mod"))
         {
           sscanf(argv[++i],"%d",&IMOD); 
           args=args+2;
         }
         else if(!strcmp(argv[i],"-split"))
         {
           sscanf(argv[++i],"%d",&N_outread);
           args=args+2;
         }
         else if(!strcmp(argv[i],"-trim"))
         {
           sscanf(argv[++i],"%d",&n_trim);
           args=args+2;
         }
         else if(!strcmp(argv[i],"-cons"))
         {
           sscanf(argv[++i],"%d",&cons_flag);
           args=args+2;
         }
      }
  
  
      num_reads = 0;
  
      Align_Process(argv,args,argc,num_reads);
      return EXIT_SUCCESS;
  
  }
  /* end of the main */
  
  /*   Subroutine to process alignment information */
  /* ====================================================  */
  void Align_Process(char **argv,int args,int argc,int nRead)
  /* ====================================================  */
  {
       int i;
  //     char **rdname;
       char *ptr,line[500],line2[500],Nameout[100],RC,motif[5] = "GATC";
       FILE *namef,*namef2,*namef3;
       long num_align=0;
       int ac,name_flag=0;
       int num_out,num_file,out_flag;
       int hit_score = 0; 
  
       RC = '+';
  
  /*   read the input fastq file         */
       num_align = 0;
       num_file  = 0;
       num_out=0;
       out_flag=1;
  
       if((namef = fopen(argv[args+1],"w")) == NULL)
       {
         printf("ERROR main:: reads group file: %s \n",argv[args]);
         exit(1);
       }
  
       if((namef3 = fopen(argv[args+2],"w")) == NULL)
       {
         printf("ERROR main:: reads group file: %s \n",argv[args]);
         exit(1);
       }
  
  /*   read the cigar line file   */
  
       if((namef2 = fopen(argv[args],"r")) == NULL)
       {
         printf("ERROR main:: alignment file 2 \n");
         exit(1);
       }
  
       num_align = 0;
       while(!feof(namef2))
       {
         int nPair=0,add_flag=0,len,hit1,hit2,m_score;
         char *pch;
        
         fgets(line,500,namef2);
         if(feof(namef2)) break;
         out_flag=1;
  //       if((num_align+1)%4==0)
  //         line[0] = '!';
  //       printf("%d %d %d %s %s",num_out,num_align,(num_align+1)%4,Nameout,line);
         if(((num_align+3)%4)==0)
         {
           memset(line2,'\0',500);
           pch = strstr(line,motif);
           if(pch != NULL)
           {
             fprintf(namef3,"HiC-reads_%010ld_1\n",num_align/4);
             fprintf(namef,"@HiC-reads_%010ld_1\n",num_align/4);
           }
           else
           {
             fprintf(namef3,"HiC-reads_%010ld_N\n",num_align/4);
             fprintf(namef,"@HiC-reads_%010ld_N\n",num_align/4);
  //           printf("%s %s\n",pch,line);
           }
           fprintf(namef,"%s",line);
         }
         else
         {
           if(((num_align)%4)!=0)
             fprintf(namef,"%s",line);
         }
         num_align++;
       }
       fclose(namef);
       fclose(namef2);
       fclose(namef3);
  
       printf("number of reads processed: %ld\n",num_align);
  
  }
  
  
  #define SWAP(a,b) temp=(a);(a)=b;(b)=temp;
  
  /*   Subroutine to sort an array arr[0,...,n-1] into ascending order while
       making the corresponding reaarangement of the array brr[0,...,n-1]
       by the use of Quicksort (Sedgwick, R. 1978, Communications o fthe ACM,
       vol. 21, pp. 847-857) also see Numerical Recipes in C                  */  
  
  /* =============================== */
  void ArraySort_Long(int n, B64_long *arr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,NSTACK=50,istack[NSTACK];
       B64_long a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
               }
               arr[i+1]=a;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /* =============================== */
  void ArraySort_Int(int n, int *arr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,NSTACK=50,istack[NSTACK];
       int a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
               }
               arr[i+1]=a;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /* =============================== */
  void ArraySort_Mix(int n, B64_long *arr, int *brr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       B64_long a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /* =============================== */
  void ArraySort_Int2(int n, int *arr, int *brr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /*   function to sort an array into a decreasing order:  a>b>c>....    */  
  /* =============================== */
  void ArraySort2_Int2(int n, int *arr, int *brr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]>=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
  
            if(arr[m]<arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
            }
  
            if(arr[m+1]<arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(arr[m]<arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            for(;;)
            {
               do i++; while (arr[i]>a);
               do j--; while (arr[j]<a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  /* =============================== */
  void ArraySort_Mix3(int n, B64_long *arr, int *brr, int *crr)
  /* =============================== */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,c,NSTACK=50,istack[NSTACK];
       B64_long a,temp,MIN=7;
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               a=arr[j];
               b=brr[j];
               c=crr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(arr[i]<=a) break;
                  arr[i+1]=arr[i];
                  brr[i+1]=brr[i];
                  crr[i+1]=crr[i];
               }
               arr[i+1]=a;
               brr[i+1]=b;
               crr[i+1]=c;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            SWAP(arr[k],arr[m+1]);
            SWAP(brr[k],brr[m+1]);
            SWAP(crr[k],crr[m+1]);
  
            if(arr[m]>arr[ir])
            {
              SWAP(arr[m],arr[ir]);
              SWAP(brr[m],brr[ir]);
              SWAP(crr[m],crr[ir]);
            }
  
            if(arr[m+1]>arr[ir])
            {
              SWAP(arr[m+1],arr[ir]);
              SWAP(brr[m+1],brr[ir]);
              SWAP(crr[m+1],crr[ir]);
            }
  
            if(arr[m]>arr[m+1])
            {
              SWAP(arr[m],arr[m+1]);
              SWAP(brr[m],brr[m+1]);
              SWAP(crr[m],crr[m+1]);
            }
  
            i=m+1;
            j=ir;
            a=arr[m+1];
            b=brr[m+1];
            c=crr[m+1];
            for(;;)
            {
               do i++; while (arr[i]<a);
               do j--; while (arr[j]>a);
               if(j<i) break;
               SWAP(arr[i],arr[j]);
               SWAP(brr[i],brr[j]);
               SWAP(crr[i],crr[j]);
            }
            arr[m+1]=arr[j];
            arr[j]=a;
            brr[m+1]=brr[j];
            brr[j]=b;
            crr[m+1]=crr[j];
            crr[j]=c;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /*   to swap the string arrays           */
  /* ============================================= */
  void s_swap(char **Pair_Name, int i, int j)
  /* ============================================= */
  {
       char temp[Max_N_NameBase];
  
       strcpy(temp,Pair_Name[i]);
       strcpy(Pair_Name[i],Pair_Name[j]);
       strcpy(Pair_Name[j],temp);
  }
  
  
  /*   to sort the string array in order          */
  /* ============================================= */
  void ArraySort_String(int n, char **Pair_Name, int *brr)
  /* ============================================= */
  {
       int i,ir=n-1,j,k,m=0,jstack=0,b,NSTACK=50,istack[NSTACK];
       int temp,MIN=7;
       char p[Max_N_NameBase];
  
       for(;;)
       {
  /*      Insertion sort when subarray is small enough    */
          if(ir-m<MIN)
          {
            for(j=m+1;j<=ir;j++)
            {
               strcpy(p,Pair_Name[j]);
               b=brr[j];
               for(i=j-1;i>=m;i--)
               {
                  if(strcmp(Pair_Name[i],p)<=0) break;
                  strcpy(Pair_Name[i+1],Pair_Name[i]);
                  brr[i+1]=brr[i];
               }
               strcpy(Pair_Name[i+1],p);
               brr[i+1]=b;
            }
            if(!jstack) return;
            ir=istack[jstack--];
            m=istack[jstack--];
          }
          else
          {
            k=(m+ir)>>1;
            s_swap(Pair_Name,k,m+1);
            SWAP(brr[k],brr[m+1]);
  
            if(strcmp(Pair_Name[m],Pair_Name[ir])>0)
            {
              s_swap(Pair_Name,m,ir);
              SWAP(brr[m],brr[ir]);
            }
  
            if(strcmp(Pair_Name[m+1],Pair_Name[ir])>0)
            {
              s_swap(Pair_Name,m+1,ir);
              SWAP(brr[m+1],brr[ir]);
            }
  
            if(strcmp(Pair_Name[m],Pair_Name[m+1])>0)
            {
              s_swap(Pair_Name,m,m+1);
              SWAP(brr[m],brr[m+1]);
            }
  
            i=m+1;
            j=ir;
            strcpy(p,Pair_Name[m+1]);
            b=brr[m+1];
            for(;;)
            {
               do i++; while (strcmp(Pair_Name[i],p)<0);
               do j--; while (strcmp(Pair_Name[j],p)>0);
               if(j<i) break;
               s_swap(Pair_Name,i,j);
               SWAP(brr[i],brr[j]);
            }
            strcpy(Pair_Name[m+1],Pair_Name[j]);
            strcpy(Pair_Name[j],p);
            brr[m+1]=brr[j];
            brr[j]=b;
            jstack+=2;
  
  /*        Push pointers to larger subarray on stack      */
  /*        process smaller subarray immediately           */
            if(jstack>NSTACK)
            {
               printf("Stack error: NSTACK too small\n");
               exit(0);
            }
            if(ir-i+1>=j-m)
            {
              istack[jstack]=ir;
              istack[jstack-1]=i;
              ir=j-1;
            }
            else
            {
              istack[jstack]=j-1;
              istack[jstack-1]=m;
              m=i;
            }
          }
       }
  }
  
  
  /* creat an int matrix with subscript ange m[nrl...nrh][ncl...nch]  */
  int     **imatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch)
  {
          B64_long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
          int  **m;
  
          /* allocate pointers to rows        */
          if((m=(int **)calloc(nrow,sizeof(int*)))==NULL)
          {
             printf("error imatrix: calloc error No. 1 \n");
             return(NULL);
          }
          m+=0;
          m-=nrl;
  
          /* allocate rows and set pointers to them        */
          if((m[nrl]=(int *)calloc(nrow*ncol,sizeof(int)))==NULL)
          {
             printf("error imatrix: calloc error No. 2 \n");
             return(NULL);
          }
          m[nrl]+=0;
          m[nrl]-=nrl;
  
          for(i=nrl+1;i<=nrh;i++)
             m[i]=m[i-1]+ncol;
          /* return pointer to array of pointers to rows   */
          return m;
  }
  
  /* creat char matrix with subscript ange cm[nrl...nrh][ncl...nch]  */
  char    **cmatrix(B64_long nrl,B64_long nrh,B64_long ncl,B64_long nch)
  {
          B64_long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
          char **cm;
  
          /* allocate pointers to rows        */
          if((cm=(char **)calloc(nrow,sizeof(char*)))==NULL)
          {
             printf("error cmatrix: calloc error No. 1 \n");
             return(NULL);
          }
          cm+=0;
          cm-=nrl;
  
          /* allocate rows and set pointers to them        */
          if((cm[nrl]=(char *)calloc(nrow*ncol,sizeof(char)))==NULL)
          {
             printf("error cmatrix: calloc error No. 2 \n");
             return(NULL);
          }
          cm[nrl]+=0;
          cm[nrl]-=nrl;
  
          for(i=nrl+1;i<=nrh;i++)
             cm[i]=cm[i-1]+ncol;
          /* return pointer to array of pointers to rows   */
          return cm;
  }
  
  