/* ========================== C MeatAxe =============================
   zqt.c - Projection on quotient.

   (C) Copyright 1993 Michael Ringe, Lehrstuhl D fuer Mathematik,
   RWTH Aachen, Germany  <mringe@tiffy.math.rwth-aachen.de>
   This program is free software; see the file COPYING for details.
   ================================================================== */


/* $Id: zquot.c,v 2.0 1993/10/14 18:54:18 mringe Exp $
 *
 * $Log: zquot.c,v $
 * Revision 2.0  1993/10/14  18:54:18  mringe
 * MeatAxe-2.0, Phase I
 *
 * Revision 2.0  1993/10/14  18:54:18  mringe
 * MeatAxe-2.0, Phase I
 *
 * Revision 1.8  1993/10/11  19:05:28  mringe
 * Neue Library-Struktur.
 *
 * Revision 1.7  1993/10/06  04:41:05  mringe
 * utils Library eliminiert.
 *
 * Revision 1.6  1993/10/05  23:35:02  mringe
 * zerrno eliminiert.
 *
 * Revision 1.5  1993/08/12  16:14:23  mringe
 * Include <string.h>.
 *
 * Revision 1.4  1993/08/06  15:09:31  mringe
 * ERR_NOECH in ERR_NOTECH umbenannt.
 *
 * Revision 1.3  1993/08/06  14:01:59  mringe
 * Neuer File-header.
 *
 * Revision 1.2  1993/02/17  11:16:12  mringe
 * Include-Files...
 *
 * Revision 1.1  1993/02/10  19:40:54  mringe
 * Initial revision
 *
 * Revision 1.1  1993/01/28  13:42:30  mringe
 * Initial revision
 *
 */

#include <string.h>
#include <stdlib.h>
#include "meataxe.h"




/* ------------------------------------------------------------------
   Global data
   ------------------------------------------------------------------ */

static long *mypiv;
static int ownpiv = 0;
static long *locn = NULL;
static long vdim, qdim, sdim;
static PTR mysubsp;


/* ------------------------------------------------------------------
   zquotinit() - Initialize. Finds the insignificant rows.

   Arguments: <subspace> must be a matrix in semi-echelon form with
	<znoc> columns and <dim> rows. <piv> must point to a pivot
	table for this matrix. If <piv>==NULL, a pivot table is
	built by a call to zmkpivot().

   Return value: Returns 0 on success, -1 on error.
   ------------------------------------------------------------------ */

int zquotinit(subspace,dim,piv)
PTR subspace;
long dim;
long *piv;

{
    long i;
    char *ispiv;

    if (locn != NULL)
    {
	free(locn);
	locn = NULL;
    }
    if (ownpiv)
    {
	free(mypiv);
	ownpiv = 0;
    }

    if (piv != NULL)
    {
    	mypiv = piv;
	ownpiv = 0;
    }
    else
    {
	mypiv = (long *)malloc((size_t)(dim+1) * sizeof (long));
	if (mypiv == NULL)
	{
	    MTXFAIL(ERR_NOMEM,-1);
	}
	ownpiv = 1;
	if (zmkpivot(subspace,dim,mypiv) != 0)
	    return -1;
    }
    vdim = znoc;
    sdim = dim;
    mysubsp = subspace;

    ispiv = (char *) malloc((size_t)(znoc+1));
    locn = (long *) malloc((size_t)(vdim-sdim+1) * sizeof(long));
    memset(ispiv,0,(size_t)(znoc+1));

    for (i = 1; i <= sdim; ++i)
    {
	if (ispiv[mypiv[i]])
	{
	    free(ispiv);
	    MTXFAIL(ERR_NOTECH,-1);
	}
	ispiv[mypiv[i]] = 1;
    }
    qdim = 0;
    for (i = 1; i <= znoc; ++i)
    {
	if (!ispiv[i])
	    locn[++qdim] = i;
    }
    free(ispiv);
    return 0;
}


/* ------------------------------------------------------------------
   zquot() - Projection of rows on the quotient.

   Parameters: <space> must be a matrix with <vdim> columns and
	<dim> rows.

   Return value: Returns a matrix with <dim> rows and <qdim>
       columns.
   ------------------------------------------------------------------ */

PTR zquot(space,dim)
PTR space;
long dim;

{
    long i, k;
    PTR q, qx, sx, tmp;

    zsetlen(zfl,qdim);
    qx = q = zalloc(dim);
    zsetlen(zfl,vdim);
    sx = space;
    tmp = zalloc((long)1);

    for (i = 1; i <= dim; ++i)
    {
	zsetlen(zfl,vdim);
	zmoverow(tmp,sx);
	zadvance(&sx,(long)1);
	zcleanrow(tmp,mysubsp,sdim,mypiv);
	zsetlen(zfl,qdim);
	zmulrow(qx,F_ZERO);
	for (k = 1; k <= qdim; ++k)
	    zinsert(qx,k,zextract(tmp,locn[k]));
	zadvance(&qx,(long)1);
    }
    free(tmp);
    zsetlen(zfl,vdim);
    return q;
}


/* ------------------------------------------------------------------
   zquotop() - Calculate the action of a square matrix on the
	quotient.

   Parameters: <matrix> must be a square matrix of dimension <vdim>.

   Return value: Returns a square matrix of dimension <qdim>.
   ------------------------------------------------------------------ */

PTR zquotop(matrix)
PTR matrix;

{
    long i, k;
    PTR q, qx, sx, tmp;

    zsetlen(zfl,qdim);
    qx = q = zalloc(qdim);
    zsetlen(zfl,vdim);
    tmp = zalloc((long)1);

    for (i = 1; i <= qdim; ++i)
    {
	zsetlen(zfl,vdim);
        sx = matrix;
	zadvance(&sx,locn[i]-1);
	zmoverow(tmp,sx);
	zcleanrow(tmp,mysubsp,sdim,mypiv);
	zsetlen(zfl,qdim);
	zmulrow(qx,F_ZERO);
	for (k = 1; k <= qdim; ++k)
	    zinsert(qx,k,zextract(tmp,locn[k]));
	zadvance(&qx,(long)1);
    }
    free(tmp);
    zsetlen(zfl,vdim);
    return q;
}