mtrans.c

IML package provides efficient routines to solve nonsingular systems of linear equations, certifie

      cblas_dgemv(CblasRowMajor, CblasNoTrans, n, n, 1.0, B, n, \
		  A+j-1, n, 0.0, C, 1);
      Dmod(p, C, 1, n, 1);
      for (i = 0; i < n-1; i++) { C[i] = fmod((p-1)*C[i], p1); }
      for (i = 0; i < n-1; i++) { for (k = 0; k < n; k++) { B[i*n+k] = 0; } }
      cblas_dger(CblasRowMajor, n-1, n, 1.0, C, 1, B+(n-1)*n, 1, B, n);
      Dmod(p, B, n-1, n, n);
      for (i = n-2; i > j-2; i--)
	{
	  ++count;
	  cblas_dswap(n, B+i*n, 1, B+(i+1)*n, 1);
	}
      if (count % 2 == 0)
	for (i = 0; i < n*n; i++) { B[i] = fmod(det*B[i], p1); }
      else
	for (i = 0; i < n*n; i++) { B[i] = fmod((p-det)*B[i], p1); }	
      { XFREE(P); XFREE(rp); XFREE(C); }
      return B;
    }
}
  



/* 
 * Calling Sequence:
 *   r/-1 <-- mBasis(p, A, n, m, basis, nullsp, B, N)
 *
 * Summary:
 *   Compute a basis for the rowspace and/or a basis for the left nullspace 
 *   of a mod p matrix
 *  
 * Description:
 *   Given a n x m mod p matrix A, the function computes a basis for the
 *   rowspace B and/or a basis for the left nullspace N of A. Row vectors in 
 *   the r x m matrix B consist of basis of A, where r is the rank of A in
 *   Z/pZ. If r is zero, then B will be NULL. Row vectors in the n-r x n
 *   matrix N consist of the left nullspace of A. N will be NULL if A is full
 *   rank.
 *
 *   The pointers are passed into argument lists to store the computed basis 
 *   and nullspace. Upon completion, the rank r will be returned. The 
 *   parameters basis and nullsp control whether to compute basis and/or
 *   nullspace. If set basis and nullsp in the way that both basis and 
 *   nullspace will not be computed, an error message will be printed and 
 *   instead of rank r, -1 will be returned.
 *  
 * Input:
 *        p: FiniteField, prime modulus
 *           if p is a composite number, the routine will still work if no 
 *           error message is returned
 *        A: 1-dim Double array length n*m, representation of a n x m mod p 
 *           matrix. The entries of A are casted from integers
 *        n: long, row dimension of A
 *        m: long, column dimension of A
 *    basis: 1/0, flag to indicate whether to compute basis for rowspace or 
 *           not
 *         - basis = 1, compute the basis
 *         - basis = 0, not compute the basis
 *   nullsp: 1/0, flag to indicate whether to compute basis for left nullspace
 *           or not
 *         - nullsp = 1, compute the nullspace
 *         - nullsp = 0, not compute the nullspace
 *
 * Output:
 *   B: pointer to (Double *), if basis = 1, *B will be a 1-dim r*m Double
 *      array, representing the r x m basis matrix. If basis = 1 and r = 0, 
 *      *B = NULL
 *  
 *   N: pointer to (Double *), if nullsp = 1, *N will be a 1-dim (n-r)*n Double
 *      array, representing the n-r x n nullspace matrix. If nullsp = 1 and
 *      r = n, *N = NULL.
 *
 * Return:
 *   - if basis and/or nullsp are set to be 1, then return the rank r of A 
 *   - if both basis and nullsp are set to be 0, then return -1
 *
 * Precondition:
 *   n*(p-1)^2 <= 2^53-1 = 9007199254740991
 *
 * Note:
 *   - In case basis = 0, nullsp = 1, A will be destroyed inplace. Otherwise,
 *     A will not be changed.
 *   - Space of B and/or N will be allocated in the function
 *
 */

long
mBasis (const FiniteField p, Double *A, const long n, const long m, \
	const long basis, const long nullsp, Double **B, Double **N)
{
  long i, r;
  long *P, *rp;
  FiniteField d[1];
  Double *A1, *B1, *N1, *U;

  P = XMALLOC(long, n+1);
  for (i = 0; i < n+1; i++) { P[i] = i; }
  rp = XCALLOC(long, n+1);
  d[0] = 1;
  if ((basis == 1) && (nullsp == 1))
    {
      A1 = XMALLOC(Double, n*m);
      DCopy(n, m, A, m, A1, m);
      RowEchelonTransform(p, A1, n, m, 1, 1, 1, 0, P, rp, d);
      r = rp[0];
      U = XCALLOC(Double, n*n);
      for (i = 0; i < n; i++) { U[i*n+i] = 1; }
      if (r != 0) { DCopy(n, r, A1, m, U, n); }
      for (i = r; i > 0; i--)
	{ if (P[i] != i) { cblas_dswap(n, U+i-1, n, U+P[i]-1, n); } }
      if (r == 0) { *B = NULL; }
      else
	{
	  *B = XMALLOC(Double, r*m);
	  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, r, m, n, \
		      1.0, U, n, A, m, 0.0, *B, m);
	  Dmod(p, *B, r, m, m);
	}
      if (r ==n) { *N = NULL; }
      else
	{
	  *N = XMALLOC(Double, (n-r)*n);
	  DCopy(n-r, n, U+r*n, n, *N, n);;
	}
      { XFREE(A1); XFREE(U); XFREE(P); XFREE(rp); }
      return r;
    }
  else if ((basis == 1) && (nullsp == 0))
    {
      A1 = XMALLOC(Double, n*m);
      DCopy(n, m, A, m, A1, m);
      RowEchelonTransform(p, A1, n, m, 1, 0, 1, 0, P, rp, d);
      r = rp[0];
      if (r == 0) 
	{
	  *B = NULL; 
	  { XFREE(A1); XFREE(P); XFREE(rp); }
	  return r;
	}
      U = XCALLOC(Double, r*n);
      DCopy(r, r, A1, m, U, n);
      for (i = r; i > 0; i--)
	{ if (P[i] != i) { cblas_dswap(r, U+i-1, n, U+P[i]-1, n); } }
      *B = XMALLOC(Double, r*m);
      cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, r, m, n, \
		  1.0, U, n, A, m, 0.0, *B, m);
      Dmod(p, *B, r, m, m);
      { XFREE(A1); XFREE(U); XFREE(P); XFREE(rp); }
      return r;
    }
  else if ((basis == 0) && (nullsp == 1))
    {
      RowEchelonTransform(p, A, n, m, 0, 1, 1, 0, P, rp, d);
      r = rp[0];
      if (r == n)
	{
	  *N == NULL;
	  { XFREE(P); XFREE(rp); }
	  return r; 
	}
      *N = XCALLOC(Double, (n-r)*n);
      if (r != 0) { DCopy(n-r, r, A+r*m, m, *N, n); }
      for (i = 0; i < n-r; i++) { (*N)[i*n+r+i] = 1; }
      for (i = r; i > 0; i--)
	{ if (P[i] != i) { cblas_dswap(n-r, *N+i-1, n, *N+P[i]-1, n); } }
      { XFREE(P); XFREE(rp); }
      return r;
    }
  else
    {
      iml_error("In mBasis, both basis and nullsp are zero.");
      return -1;
    }
}



/*
 * Calling Sequence:
 *   det <-- mDeterminant(p, A, n)
 * 
 * Summary:
 *   Compute the determinant of a square mod p matrix
 *
 * Input:
 *   p: FiniteField, prime modulus
 *      if p is a composite number, the routine will still work if no error 
 *      message is returned
 *   A: 1-dim Double array length n*n, representation of a n x n mod p matrix.
 *      The entries of A are casted from integers
 *   n: long, dimension of A
 *
 * Output:
 *   det(A) mod p, the determinant of square matrix A
 *
 * Precondition:
 *   ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
 *
 * Note:
 *   A is destroyed inplace
 *
 */
long 
mDeterminant (const FiniteField p, Double *A, const long n)
{
  long i, count=0;
  long *P, *rp;
  FiniteField det, d[1];

  P = XMALLOC(long, n+1);
  for (i = 0; i < n+1; i++) { P[i] = i; }
  rp = XCALLOC(long, n+1);
  d[0] = 1;
  RowEchelonTransform(p, A, n, n, 0, 0, 0, 1, P, rp, d);
  det = d[0];
  if (det != 0)
    {
      for (i = 1; i < n+1; i++)	{ if (P[i] != i) { ++count; } }
      if (count % 2 == 0) 
	{ 
	  { XFREE(P); XFREE(rp); }
	  return det;
	}
      else
	{
	  { XFREE(P); XFREE(rp); }
	  return p-det;
	}
    }
  { XFREE(P); XFREE(rp); }     
  return det;
}


/*
 * Calling Sequence:
 *   1/0 <-- mInverse(p, A, n)
 * 
 * Summary:
 *   Certified compute the inverse of a mod p matrix inplace
 *
 * Description:
 *   Given a n x n mod p matrix A, the function computes A^(-1) mod p 
 *   inplace in case A is a nonsingular matrix in Z/Zp. If the inverse does
 *   not exist, the function returns 0.
 *
 *   A will be destroyed at the end in both cases. If the inverse exists, A is
 *   inplaced by its inverse. Otherwise, the inplaced A is not the inverse.
 *
 * Input:
 *   p: FiniteField, prime modulus
 *      if p is a composite number, the routine will still work if no error 
 *      message is returned
 *   A: 1-dim Double array length n*n, representation of a n x n mod p matrix.
 *      The entries of A are casted from integers
 *   n: long, dimension of A
 *
 * Return: 
 *   - 1, if A^(-1) mod p exists
 *   - 0, if A^(-1) mod p does not exist
 *
 * Precondition:
 *   ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
 *
 * Note:
 *   A is destroyed inplace
 *
 */

long
mInverse (const FiniteField p, Double *A, const long n)
{
  long i;
  long *P, *rp;
  FiniteField d[1];

  P = XMALLOC(long, n+1);
  for (i = 0; i < n+1; i++) { P[i] = i; }
  rp = XCALLOC(long, n+1);
  d[0] = 1;
  RowEchelonTransform(p, A, n, n, 1, 1, 1, 0, P, rp, d);
  if (rp[0] == n) 
    {
      for (i = n; i > 0; i--)
	if (P[i] != i)
	  cblas_dswap(n, A+i-1, n, A+P[i]-1, n);
      { XFREE(P); XFREE(rp); }
      return 1;
    }
  { XFREE(P);  XFREE(rp); }
  return 0;
}




/*
 * Calling Sequence:
 *   r <-- mRank(p, A, n, m)
 *
 * Summary:
 *   Compute the rank of a mod p matrix
 *
 * Input:
 *   p: FiniteField, prime modulus
 *      if p is a composite number, the routine will still work if no 
 *      error message is returned
 *   A: 1-dim Double array length n*m, representation of a n x m mod p 
 *      matrix. The entries of A are casted from integers
 *   n: long, row dimension of A
 *   m: long, column dimension of A
 *   
 * Return:
 *   r: long, rank of matrix A
 *
 * Precondition:
 *   ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
 *
 * Note:
 *   A is destroyed inplace
 *
 */

long 
mRank (const FiniteField p, Double *A, const long n, const long m)
{
  long i, r;
  long *P, *rp;
  FiniteField d[1];

  P = XMALLOC(long, n+1);
  for (i = 0; i < n+1; i++) { P[i] = i; }
  rp = XCALLOC(long, n+1);
  d[0] = 1;
  RowEchelonTransform(p, A, n, m, 0, 0, 0, 0, P, rp, d);
  r = rp[0];
  { XFREE(P); XFREE(rp); }
  return r;
}




/*
 * Calling Sequence:
 *   rp <-- mRankProfile(p, A, n, m)
 *
 * Summary:
 *   Compute the rank profile of a mod p matrix
 *
 * Input:
 *   p: FiniteField, prime modulus
 *      if p is a composite number, the routine will still work if no 
 *      error message is returned
 *   A: 1-dim Double array length n*m, representation of a n x m mod p 
 *      matrix. The entries of A are casted from integers
 *   n: long, row dimension of A
 *   m: long, column dimension of A
 *   
 * Return:
 *   rp: 1-dim long array length n+1, where
 *     - rp[0] is the rank of matrix A
 *     - rp[1..r] is the rank profile of matrix A
 *
 * Precondition:
 *   ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
 *
 * Note:
 *   A is destroyed inplace
 *
 */

long *
mRankProfile (const FiniteField p, Double *A, const long n, const long m)
{
  long i;
  long *P, *rp;
  FiniteField d[1];

  P = XMALLOC(long, n+1);
  for (i = 0; i < n+1; i++) { P[i] = i; }
  rp = XCALLOC(long, n+1);
  d[0] = 1;
  RowEchelonTransform(p, A, n, m, 0, 0, 0, 0, P, rp, d);
  XFREE(P);
  return rp;
}