/******************************************************************************
 *
 *       ELMER, A Computational Fluid Dynamics Program.
 *
 *       Copyright 1st April 1995 - , Center for Scientific Computing,
 *                                    Finland.
 *
 *       All rights reserved. No part of this program may be used,
 *       reproduced or transmitted in any form or by any means
 *       without the written permission of CSC.
 *
 ******************************************************************************/

/*******************************************************************************
 *
 *     MATC operator routines.
 *
 *******************************************************************************
 *
 *                     Author:       Juha Ruokolainen
 *
 *                    Address: Center for Scientific Computing
 *                                Tietotie 6, P.O. BOX 405
 *                                  02101 Espoo, Finland
 *                                  Tel. +358 0 457 2723
 *                                Telefax: +358 0 457 2302
 *                              EMail: Juha.Ruokolainen@csc.fi
 *
 *                       Date: 30 May 1996
 *
 *                Modified by:
 *
 *       Date of modification:
 *
 ******************************************************************************/
/***********************************************************************
|
|   Oper.C - Last Edited 15. 8. 1988
|
***********************************************************************/

/*======================================================================
|Syntax of the manual pages:
|
|FUNCTION NAME(...) params ...
|
$  usage of the function and type of the parameters
?  explane the effects of the function
=  return value and the type of value if not of type int
@  globals effected directly by this routine
!  current known bugs or limitations
&  functions called by this function
~  these functions may intrest you as an alternative function or
|  because they control this function somehow
=====================================================================*/


/*
 * $Id: oper.c,v 1.1.1.1 2005/04/14 13:29:14 vierinen Exp $ 
 *
 * $Log: oper.c,v $
 * Revision 1.1.1.1  2005/04/14 13:29:14  vierinen
 * initial matc automake package
 *
 * Revision 1.2  1998/08/01 12:34:51  jpr
 *
 * Added Id, started Log.
 * 
 *
 */

#include "elmer/matc.h"

#ifdef TYPE
#undef TYPE
#endif
#ifdef NROW
#undef NROW
#endif
#ifdef NCOL
#undef NCOL
#endif
#ifdef MATR
#undef MATR
#endif
#ifdef M
#undef M
#endif
#ifdef MATSIZE
#undef MATSIZE
#endif

#define TYPE(mat) (mat)->type
#define NROW(mat) (mat)->nrow
#define NCOL(mat) (mat)->ncol
#define MATR(mat) (mat)->data
#define MATSIZE(mat) (NROW(mat) * NCOL(mat) * sizeof(double))

#define MA(i,j) a[ncola * (i) + (j)]
#define MB(i,j) b[ncolb * (i) + (j)]
#define MC(i,j) c[ncolc * (i) + (j)]

MATRIX *mat_new(type, nrow, ncol) int type, nrow, ncol;
{
   MATRIX *res;

    res = (MATRIX *)ALLOCMEM(MATRIXSIZE);
    TYPE(res) = type;
    NROW(res) = nrow;
    NCOL(res) = ncol;
    MATR(res) = (double *)ALLOCMEM(MATSIZE(res));

    return res;
}

MATRIX *mat_copy(mat) MATRIX *mat;
{
    MATRIX *res;

    if (mat == (MATRIX *)NULL) return NULL;
     
    res = mat_new(TYPE(mat), NROW(mat), NCOL(mat));
    memcpy((char *)MATR(res), (char *)MATR(mat), MATSIZE(mat));

    return res;
}

void mat_free(mat) MATRIX *mat;
{
    if (mat == (MATRIX *)NULL) return;

    FREEMEM((char *)MATR(mat));
    FREEMEM((char *)mat);
}

MATRIX *opr_vector(A, B)
       MATRIX *A, *B;
{
  MATRIX *C;

  int i, n, inc;

  double *a = MATR(A), *b = MATR(B), *c;

  inc = ((int)*b > (int)*a) ? 1:(-1);
  n = abs((int)*b - (int)*a) + 1;

  C = mat_new(TYPE_DOUBLE, 1, n);
  c = MATR(C);

  for(i = 0; i < n; i++) 
    *c++ = (int)*a + i*inc;

  return C;
} 

MATRIX *opr_resize(A, B)
       MATRIX *A, *B;
{
  MATRIX *C;

  double *a = MATR(A), *b = MATR(B), *c;
  int i, j, n, m;
  
  if (NCOL(B) >= 2)
  {
    i = *b++; j = *b++;
  } 
  else
  {
    i = 1; j = *b;
  }
  
  if (i < 1 || j < 1)
    error("resize: invalid size for and array");
  
  C = mat_new(TYPE(A), i, j);
  c = MATR(C);

  n = i * j; m = NROW(A) * NCOL(A);

  for(i = j = 0; i < n; i++)
  {
    *c++ = a[j++];
    if (j == m) j = 0;
  }

  return C;
}

MATRIX *opr_apply(A)
     MATRIX *A;
{
  VARIABLE *store, *ptr;
  MATRIX *C = NULL;
  
  store = var_temp_new(TYPE_STRING, NROW(A), NCOL(A));
  REFCNT(store) = 0;
  mat_free(store->this);
  store->this = A;
  REFCNT(store)++;
  
  ptr = (VARIABLE *)com_apply(store);

  var_delete_temp(store);

  if ( ptr ) C = mat_copy(ptr->this);

  return C;
}

MATRIX *opr_add(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  int i;

  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  double value;
  
  if (nrowa == nrowb && ncola == ncolb)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C); 
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++)  *c++ = *a++ + *b++;
  }

  else if (nrowa == 1 && ncola == 1)
  {
    C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
    value = *a; nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = value + *b++;
  }
  
  else if (nrowb == 1 && ncolb == 1)
  {
    C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
    value = *b; nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = value + *a++;
  }

  else
    error("Add: Incompatible for addition.\n");
  
  return C;
}

MATRIX *opr_minus(A)
     MATRIX *A;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);

  double *a = MATR(A), *c;

  C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);

  nrowa *= ncola;
  for(i = 0; i < nrowa; i++) *c++ = -*a++;
  
  return C;
}

MATRIX *opr_subs(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  double value;

  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;
  
  if (nrowa == nrowb && ncola == ncolb)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = *a++ - *b++;
  }
  else if (nrowa == 1 && ncola == 1)
  {
    C = mat_new(TYPE(B),nrowb, ncolb); c = MATR(C);
    value = *a; nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = value - *b++;
  }
  else if (nrowb == 1 && ncolb == 1)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    value = *b; nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = *a++ - value;
  }
  else
    error("Substr: Incompatible for addition.\n");
  
  return C;
}

MATRIX *opr_mul(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  double value,s;
  int i, j, k;

  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if (nrowa == 1 && ncola == 1)
  {
    C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
    value = *a; nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = value * *b++;
  }
  else if (nrowb == 1 && ncolb == 1)
  {
    C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
    value = *b; nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = value * *a++;
  }
  else if (ncola == nrowb)
  {
    C = mat_new(TYPE(A), nrowa,ncolb); 
    c = MATR(C);
    for(i = 0; i < nrowa; i++)
    {
      for(j = 0; j < ncolb; j++)
      {
        s = 0;
	for(k = 0; k < ncola; k++) s += a[k] * MB(k,j);
	*c++ = s;
      }
      a += ncola;
    }
  }
  else if ( ncola == ncolb && nrowa == nrowb )
  {
    C = mat_new(TYPE(A), nrowa,ncolb); 
    c = MATR(C);
    k = 0;
    for( i = 0; i < nrowa; i++ )
        for( j = 0; j < ncolb; j++,k++ ) c[k] = a[k] * b[k];
  }
  else 
  {
     error("Mul: Incompatible for multiplication.\n");
  }
 
  return C;
}

MATRIX *opr_pmul(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  double value;
  int i;

  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);
  
  double *a = MATR(A), *b = MATR(B), *c;

  if (nrowa == nrowb && ncola == ncolb)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = *a++ * *b++;
  }
  else if (nrowa == 1 && ncola == 1)
  {
    C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
    value = *a; nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = value * *b++;
  }
  else if (nrowb == 1 && ncolb == 1)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    value = *b; nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = *a++ * value;
  }
  else 
    error("PMul: Incompatible for pointwise multiplication.\n");
  
  return C;
}

MATRIX *opr_div(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  double value;
  int i;

  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);
  
  double *a = MATR(A), *b = MATR(B), *c;

  if (nrowa == nrowb && ncola == ncolb)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = *a++ / *b++;
  }
  else if (nrowa == 1 && ncola == 1)
  {
    C = mat_new(TYPE(B), nrowb, ncolb); c = MATR(C);
    value = *a; nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = value / *b++;
  }
  else if (nrowb == 1 && ncolb == 1)
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    value = *b; nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = *a++ / value;
  }
  else 
    error("Div: Incompatible for division.\n");
  
  return C;
}

MATRIX *opr_pow(A,B)
     MATRIX *A, *B;
{
  MATRIX *C;

  int i, j, k, l, power;

  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);
  int ncolc;

  double *a = MATR(A), *b = MATR(B), *c;

  double *v, value;

  if (nrowb != 1 || ncolb != 1) error("Pow: Matrix ^ Matrix ?.\n");

  if (nrowa == 1 || ncola != nrowa)
  {
    C = mat_new(TYPE(A), nrowa, ncola); c = MATR(C);
    value = *b; nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = pow(*a++, value);

    return C;
  }

  power = (int)*b;

  if (power == 0) {
    C = mat_new(TYPE(A), nrowa, ncola);
    ncolc = ncola; c = MATR(C);
    for(i = 0; i < nrowa; i++) MC(i,i) = 1;
  }
  else if (abs(power) == 1)
  {
    C = mat_copy(A);
  }
  else
  {
    v = (double *)ALLOCMEM(nrowa * sizeof(double));

    C = mat_new(TYPE(A), nrowa, nrowa);
    c = MATR(C); b = MATR(A);

    for(l = 1; l < abs(power); l++) {
      for(i = 0; i < nrowa; i++)
      {
	for(j = 0; j < nrowa; j++)
	{
	  v[j] = 0.0;
	  for(k = 0; k < nrowa; k++) v[j] +=  b[k] * MA(k,j);
	}
	for(j = 0; j < nrowa; j++) *c++ = v[j];
	b += nrowa;
      }
      a = MATR(A); b = c = MATR(C);
    }
    FREEMEM((char *)v);
  }

  if (power < 0)
  {
    VARIABLE *inv, *tmp;
    tmp = (VARIABLE *)ALLOCMEM(VARIABLESIZE);
    tmp -> this = C;
    inv = mtr_inv(tmp);
    mat_free(C);
    FREEMEM((char *)tmp);
    C = inv->this;
    REFCNT(inv)++;
    var_delete_temp(inv);
  }

  return C;
}

MATRIX *opr_trans(A)
     MATRIX *A;
{
  MATRIX *C;

  int i, j;

  int ncola = NCOL(A), nrowa = NROW(A);
  int ncolc;  

  double *a = MATR(A), *c;

  C = mat_new(TYPE(A),ncola,nrowa);
  c = MATR(C); ncolc = nrowa;

  for(i = 0; i < nrowa; i++)
    for(j = 0; j < ncola; j++) MC(j,i) = *a++;
  
  return C;
}

MATRIX *opr_reduction(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  int i;

  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++, a++) *c++ = (*b++) ? (*a) : (0);
  }
  else
    error("Incompatible for reduction.\n");

  return C;
}


MATRIX *opr_lt(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb, ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++,c++) if (*a < b[i]) *c = 1;
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] < *b) *c = 1;
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] < b[i]) *c = 1;
  }
  else
    error("lt: Incompatible for comparison.\n");
  
  return C;
}


MATRIX *opr_le(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;

  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++,c++) if (*a <= b[i]) *c = 1;
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] <= *b) *c = 1;
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] <= b[i]) *c = 1;
  }
  else
    error("le: Incompatible for comparison.\n");
  
  return C;
}


MATRIX *opr_gt(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++,c++) if (*a > b[i]) *c = 1;
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] > *b) *c = 1;
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] > b[i]) *c = 1;
  }
  
  else
    error("gt: Incompatible for comparison.\n");
  
  return C;
}


MATRIX *opr_ge(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowa; i++,c++) if (*a >= b[i]) *c = 1;
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] >= *b) *c = 1;
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] >= b[i]) *c = 1;
  }
  
  else
    error("ge: Incompatible for comparison.\n");
  
  return C;
}


MATRIX *opr_eq(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++,c++) if (*a == b[i]) *c = 1;
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] == *b) *c = 1;
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] == b[i]) *c = 1;
  }
  
  else 
    error("eq: Incompatible for comparison.\n");
  
  return C;
}


MATRIX *opr_neq(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++,c++) if (*a != b[i]) *c = 1;
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] != *b) *c = 1;
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++,c++) if (a[i] != b[i]) *c = 1;
  }
  
  else 
    error("neq: Incompatible for comparison.\n");
  
  return C;
}

MATRIX *opr_or(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = (*a) || (b[i]);
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = (a[i]) || (*b);
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = (a[i]) || (b[i]);
  }
  
  else 
    error("or: Incompatible for comparison.\n");
  
  return C;
}

MATRIX *opr_and(A, B)
     MATRIX *A, *B;
{
  MATRIX *C;
  int i;
  
  int nrowa = NROW(A), ncola = NCOL(A);
  int nrowb = NROW(B), ncolb = NCOL(B);

  double *a = MATR(A), *b = MATR(B), *c;

  if ((nrowa == 1) && (ncola == 1))
  {
    C = mat_new(TYPE(B),nrowb,ncolb); c = MATR(C);
    nrowb *= ncolb;
    for(i = 0; i < nrowb; i++) *c++ = (*a) && (b[i]);
  }
  
  else if ((nrowb == 1) && (ncolb == 1))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = (a[i]) && (*b);
  }
  
  else if ((nrowa == nrowb) && (ncola == ncolb))
  {
    C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
    nrowa *= ncola;
    for(i = 0; i < nrowa; i++) *c++ = (a[i]) && (b[i]);
  }
  
  else 
    error("and: Incompatible for comparison.\n");
  
  return C;
}

MATRIX *opr_not(A)
     MATRIX *A;
{
  MATRIX *C;
  int i;

  int nrowa = NROW(A), ncola = NCOL(A);

  double *a = MATR(A), *c;

  C = mat_new(TYPE(A),nrowa,ncola); c = MATR(C);
  nrowa *= ncola;
  for(i = 0; i < nrowa; i++,c++) if (*a == 0) *c = 1;

  return C;
}