/* ========================== C MeatAxe ============================= zsp.c - Spin up and split. (C) Copyright 1993 Michael Ringe, Lehrstuhl D fuer Mathematik, RWTH Aachen, Germany This program is free software; see the file COPYING for details. ================================================================== */ /* $Id: zsp.c,v 2.6 1994/03/13 13:27:01 mringe Exp $ * * $Log: zsp.c,v $ * Revision 2.6 1994/03/13 13:27:01 mringe * Maschinenunabhaengiges Format fuer Permutationen. * * Revision 2.6 1994/03/13 13:27:01 mringe * Maschinenunabhaengiges Format fuer Permutationen. * * Revision 2.5 1994/02/13 18:26:56 mringe * Neu: os.c, os.h. * * Revision 2.4 1993/12/08 11:33:02 mringe * Neue CPU time - Funktionen. * * Revision 2.3 1993/10/21 21:57:35 mringe * Permutationen. * * Revision 2.2 1993/10/21 08:37:53 mringe * Compiler warnings. * * Revision 2.1 1993/10/20 18:17:07 mringe * MeatAxe-2.0, Phase II. * * Revision 2.0 1993/10/14 18:54:18 mringe * MeatAxe-2.0, Phase I * * Revision 1.39 1993/10/11 19:05:28 mringe * Neue Library-Struktur. * * Revision 1.38 1993/10/06 04:41:05 mringe * utils Library eliminiert. * * Revision 1.37 1993/08/24 12:50:16 mringe * Kleine Korrektur am helptext. * * Revision 1.36 1993/08/24 12:48:34 mringe * Option -T * * Revision 1.35 1993/08/24 08:40:43 mringe * Bilde Filenamen mit '.', wie CHOP. MAXDIM entfernt. * * Revision 1.34 1993/08/10 14:51:42 mringe * Include string.h * * Revision 1.33 1993/08/06 14:01:59 mringe * Neuer File-header. * * Revision 1.32 1993/08/05 16:14:53 mringe * Prototype von err() korrigiert. * * Revision 1.31 1993/08/05 15:48:54 mringe * Neues message.c * * Revision 1.30 1993/07/28 13:34:49 mringe * *** empty log message *** * * Revision 1.29 1993/07/23 13:46:27 mringe * OS-Symbole neu (SYS_xxx) * * Revision 1.28 1993/07/20 06:18:06 mringe * GAP-Output: benutze MeatAxe-Record. * * Revision 1.27 1993/07/17 19:13:05 mringe * Aenderungen fuer Borland C. * * Revision 1.26 1993/07/17 08:38:14 mringe * msg_t durch message_t ersetzt * * Revision 1.25 1993/07/13 20:30:59 mringe * Neue File i/o library. * * Revision 1.24 1993/07/08 10:09:42 mringe * Belege nur so viel Speicher, wie mit -d angegeben. * * Revision 1.23 1993/07/02 10:49:48 mringe * Deklaration von err() fuer Nicht-ANSI-Compiler * * Revision 1.22 1993/04/07 08:21:56 mringe * Verbesserte Fehlermeldung bei inkompatiblen Files * * Revision 1.21 1993/02/17 11:16:12 mringe * Include-Files... * * Revision 1.20 1993/02/10 19:40:54 mringe * Libraries angelegt (YYY und ZZZ). * * Revision 1.19 1993/01/21 10:59:00 mringe * fflush(stdout) bei verbose output. * * Revision 1.18 1993/01/11 16:08:37 hiss * Pruefe, ob seed vector und Permutationen kompatibel sind. * * Revision 1.17 1993/01/08 12:57:59 hiss * Berechne Operation auf dem Teilraum bei Permutationen. * * Revision 1.16 1993/01/08 12:31:46 hiss * Option -v (verbose). * * Revision 1.15 1993/01/06 20:59:57 mringe * getopt in zgetopt() umbenannt. * * Revision 1.14 1992/11/04 09:45:31 mringe * Behandlung von Permutationen verbessert. * * Revision 1.13 1992/11/04 09:11:47 mringe * Erlaube Permutationen als Erzeuger. * * Revision 1.12 1992/10/02 10:05:53 mringe * Benutze pseed.c zum Generieren der Vektoren. * * Revision 1.11 1992/09/27 08:58:11 mringe * help() verbessert. * * Revision 1.10 1992/09/27 08:50:57 mringe * help() ge"andert. * * Revision 1.9 1992/08/14 10:52:03 mringe * Bug bei -ss korrigiert. * * Revision 1.8 1992/07/15 13:32:01 mringe * Prototypes mit _PLi, benutze args.c * * Revision 1.7 1992/07/15 13:02:09 mringe * Optionen -sg und -sr * * Revision 1.6 1992/07/15 09:25:55 mringe * Some minor changes. * * Revision 1.5 1992/06/29 16:55:06 mringe * Hilfstext verk"urzt. * Fehler bei der Auswertung der -s Option behoben. * * Revision 1.4 1992/06/28 10:40:34 mringe * Syntax vereinfacht. * * Revision 1.3 1992/06/28 08:59:12 mringe * Neue Version. * * interactive() eingebaut * * Revision 1.1 1992/05/26 18:33:33 mringe * Initial revision * */ #include #include #include #include "meataxe.h" /* ------------------------------------------------------------------ Function prototypes ------------------------------------------------------------------ */ void parseargs _PL((int argc, char *argv[])); int main _PL((int argc, char *argv[])); void init _PL((void)); void err _PL((int code)); void try_split _PL((void)); int spin_up _PL((void)); void do_split _PL((void)); void wrmod _PL((long cnt)); /* ------------------------------------------------------------------ Constants ------------------------------------------------------------------ */ #define EXIT_SPLIT 1 /* Return value if space was split */ #define EXIT_NOSPLIT 0 /* Return value if space was not split */ #define MAXGEN 10 /* Maximal number of generators */ /* ------------------------------------------------------------------ Global data ------------------------------------------------------------------ */ enum {DONTKNOW,MATRIX,PERMUTATION} gentype = DONTKNOW; /* Type of generators */ enum {GENERATE,READ,READALL} smode = GENERATE; /* Seed mode */ int opt_G = 0; /* GAP output */ int split = 1; /* Split after spin-up */ int output = 1; /* Enable file ouput */ long seedvecno = 0; /* Selected vector (0 = all) */ long iseed = 0; /* Current vector */ long maxdim = -1; /* Limit on susbspace dimension */ int nmat = 2; /* Number of matrices to split */ int ngen = 2; /* Number of generators */ char *genname[MAXGEN] = {"G1","G2"}; char *subname[MAXGEN] = {"P4","P5"}; char *quotname[MAXGEN] = {"P6","P7"}; char *seedname = "G3"; char *vname = "V"; FILE *seedfile; PTR gen[MAXGEN]; /* Matrices to split */ PTR m3, /* Subspace */ tmp1, tmp2; /* Two extra rows */ long nor, /* Degree of representation */ dim, /* Dimension of invariant subspace */ ssdim, /* Dimension of seed vector space */ fl; /* Field order */ long *piv; /* Pivot table for the subspace */ char *ins; /* Used to find the insignificant columns */ /* ------------------------------------------------------------------ help() ------------------------------------------------------------------ */ static char *helptext[] = { "SYNTAX", " zsp [] [ [{ | }]]", " zsp [] -g [.] [ [ [{_{| }]]]",
"",
"OPTIONS",
" -s Select seed for spin-up:",
" -sg[] Generate seed vectors (Default)",
" -sr[] Read seed vectors one by one",
" -ss Read a whole space and make its closure",
" , if specified, selects a single vector",
" -n Spin up and write out invariant subspace, but do not split",
" -o No output to files, messages only",
" -V Verbose",
" -Q Quiet, no messages",
" -d Set limit on dimension of subspace",
" -T Set limit on CPU time (in seconds)",
" -g [.] Read matrices and treat the first as generators",
" (default: =). Read matrices from .i (i=1..),",
" and write output to _{.i and .i.",
"",
"FILES",
" i Generators: square matrices of equal size",
" i Seed vectors: a matrix with NOC(Seed)=NOC(G1)",
" o Action on subspace",
" o Action on quotient",
" o Invariant subspace (with -n)",
NULL};

static proginfo_t pinfo =
{ "zsp", "Spin Up And Split", "$Revision: 2.6 $", helptext };

/* ------------------------------------------------------------------
parseargs()
------------------------------------------------------------------ */

static void setname(name,i,template)
char **name;
int i;
char *template;

{
char *c;

c = malloc(strlen(template) + 10);
sprintf(c,"%s.%d",template,i+1);
name[i] = c;
}

void parseargs(argc, argv)
int argc;
char *argv[];

{ int k, j;
int opt_g = 0;

/* Parse command line
------------------ */
initargs(argc, argv, &pinfo);
while ((k = zgetopt("ionvGd:g:s:")) != OPT_END)
{ switch (k)
{
case 'o':
output = 0;
break;
case 'G':
msg_level = -100;
opt_G = 1;
break;
case 'n':
split = 0;
break;
case 'd':
maxdim = getint();
if (maxdim <= 0 || *opt_text_ptr != 0)
errexit(ERR_OPTION,"-d");
break;
case 'g':
opt_g = 1;
nmat = (int) getint();
if (*opt_text_ptr == '.')
{ ++opt_text_ptr;
ngen = (int) getint();
}
else
ngen = nmat;
if (ngen > nmat || ngen < 2 || nmat > MAXGEN ||
*opt_text_ptr != 0)
errexit(ERR_OPTION,"-g");
break;
case 's': /* Seed mode */
switch (opt_text[0])
{ case 'g': smode = GENERATE; ++opt_text_ptr; break;
case 'r': smode = READ; ++opt_text_ptr; break;
case 's': smode = READALL; ++opt_text_ptr; break;
default: break; /* Allow -s for -sg */
}
if (smode == GENERATE || smode == READ)
{ if (*opt_text_ptr != 0)
{ seedvecno = getint();
if (seedvecno <= 0) errexit(ERR_OPTION,"-s");
}
}
if (*opt_text_ptr != 0) errexit(ERR_OPTION,"-s");
break;
}
}

/* Set file names
-------------- */
k = opt_ind; /* Index of first file arg */
if (!opt_g) /* Old Syntax */
{ switch (argc - k)
{ case 0:
break;
case 7:
quotname[0] = argv[k+5];
quotname[1] = argv[k+6];
case 5:
subname[0] = argv[k+3];
subname[1] = argv[k+4];
if (!split) errexit(ERR_BADUSAGE,"zsp");
case 4:
if (!split)
vname = argv[k+3];
case 3:
genname[0] = argv[k+0];
genname[1] = argv[k+1];
seedname = argv[k+2];
break;
default:
errexit(ERR_BADUSAGE,"zsp");
}
}
else /* The -g option was used */
{
for (j = 0; j < nmat; ++j)
{
setname(genname,j,"G");
setname(subname,j,"S");
setname(quotname,j,"Q");
}
seedname = "G0";
switch (argc - k)
{ case 4:
if (!split) errexit(ERR_BADUSAGE,"zsp");
for (j = 0; j < nmat; ++j)
setname(quotname,j,argv[k+3]);
case 3:
if (!split)
vname = argv[k+2];
else
for (j = 0; j < nmat; ++j)
setname(subname,j,argv[k+2]);
case 2:
seedname = argv[k+1];
case 1:
for (j = 0; j < nmat; ++j)
setname(genname,j,argv[k]);
case 0:
break;
default:
errexit(ERR_BADUSAGE,"zsp");
}
}
}

/* ------------------------------------------------------------------
main()
------------------------------------------------------------------ */

int main(argc, argv)
int argc;
char *argv[];

{
PTR base;

mtxinit();
parseargs(argc,argv);
init();

switch (smode)
{
case GENERATE:
base = zalloc(ssdim);
zreadvec(seedfile,base,ssdim);
zpseed_init(ssdim,base);
if (seedvecno == 0)
{
iseed = zpseed_make((long)1);
while (iseed >= 0)
{
zmoverow(m3,zpseed_vec);
try_split();
iseed = zpseed_next();
}
}
else
{
iseed = zpseed_make(seedvecno);
if (iseed < 0) err('v');
zmoverow(m3,zpseed_vec);
try_split();
}
break;
case READ:
if (seedvecno == 0)
{
for (iseed = 1; iseed <= ssdim; ++iseed)
{
zreadvec(seedfile,m3,(long)1);
try_split();
}
}
else
{ if (seedvecno > ssdim) err('v');
zseek(seedfile,seedvecno);
zreadvec(seedfile,m3,(long)1);
iseed = seedvecno;
try_split();
}
break;

case READALL:
zreadvec(seedfile,m3,ssdim);
try_split();
break;
}

/* The representation was not split. Print an
appropriate message depending on the options
used and exit with the right exit code.
--------------------------------------------- */
if (!split)
return EXIT_OK;
if (seedvecno == 0 && smode != READALL)
MESSAGE(0,("EVERY SEED VECTOR "));
MESSAGE(0,("SPANS WHOLE SPACE\n"));
return EXIT_NOSPLIT;
}

/* ------------------------------------------------------------------
init() - Read input files and initialize
------------------------------------------------------------------ */

static char wm1[] =
"zsp: Warning: Cannot calculate quotient for Permutations\n";
static char wm2[] =
"zsp: Warning: Only one permutation read from %s\n";

void init()

{
long fl2, noc, nor2, noc2;
int j;
FILE *f;

/* Open the file
-------------------- */
if (MSG2) printf("Open seed file (%s)\n",seedname);
if ((seedfile = zreadhdr(seedname,&fl,&ssdim,&noc)) == NULL)
errexit(-1,seedname);
if (MSG1) printf("Found %ld seed vectors\n",ssdim);
if (fl < 2) errexit(ERR_NOTMATRIX,seedname);
nor = noc;

/* Read the generators
------------------- */
for (j = 0; j < nmat; ++j) /* Read the generators */
{
if ((f = zreadhdr(genname[j],&fl2,&nor2,&noc2)) == NULL)
errexit(-1,genname[j]);
if (fl2 == -1) /* Permutations */
{
if (gentype == MATRIX)
errexit(ERR_NOTPERM,genname[j]);
if (gentype == DONTKNOW && split) fprintf(stderr,wm1);
gentype = PERMUTATION;
if (noc2 > 1)
fprintf(stderr,wm2,genname[j]);
if (nor2 != noc)
{
char buf[200];
sprintf(buf,"%s and %s",genname[j],seedname);
errexit(ERR_INCOMPAT,buf);
}
gen[j] = (PTR) malloc(sizeof(long)*noc);
if (gen[j] == NULL) errexit(ERR_NOMEM,"zsp");
if (zreadlong(f,(long*) gen[j],noc) != noc)
errexit(-1,genname[j]);
}
else
{
if (gentype == PERMUTATION)
errexit(ERR_NOTMATRIX,genname[j]);
gentype = MATRIX;
if (fl2 != fl || noc2 != noc)
{
char buf[200];
sprintf(buf,"%s and %s",genname[j],seedname);
errexit(ERR_INCOMPAT,buf);
}
if (noc2 != nor2) errexit(ERR_NOTSQUARE,genname[j]);
zsetlen(fl,nor);
gen[j] = zalloc(nor);
zreadvec(f,gen[j],nor);
}
fclose(f);
}

/* Set default value of maxdim
--------------------------- */
if (maxdim == -1) maxdim = nor-1;
if (smode == READALL && maxdim < ssdim) err('S');

/* Allocate memory
--------------- */
zsetlen(fl,nor);
m3 = zalloc(maxdim + 2);
tmp1 = zalloc((long)1);
tmp2 = zalloc((long)1);
piv = (long *)malloc((size_t)(nor+1)*sizeof(long));
ins = (char *)malloc((size_t)(nor+1));
if (piv == NULL || ins == NULL)
errexit(ERR_NOMEM,"zsp");
}

/* ------------------------------------------------------------------
err() - Terminate program with error message
------------------------------------------------------------------ */

void err(code)
int code;

{
fprintf(stderr,"ZSP ERROR - ");
switch (code)
{
case 'b':
fprintf(stderr,"INTERNAL ERROR\n");
break;
case '0':
fprintf(stderr,"CANNOT SPIN UP ZERO VECTOR\n");
break;
case 'v':
fprintf(stderr,"CANNOT MAKE THIS SEED VECTOR\n");
break;
case 'S':
fprintf(stderr,
"SEED SPACE IS GREATER THAN DIMENSION LIMIT");
break;
}
exit(EXIT_ERR);
}

/* ------------------------------------------------------------------
try_split()
------------------------------------------------------------------ */

void try_split()

{
int hassplit;

hassplit = spin_up();
if (opt_G)
printf("MeatAxe.SplitDimensions := [%ld,%ld];\n",dim,nor-dim);

if (!split)
{
if (smode != READALL)
MESSAGE(0,("VECTOR %5ld: ",iseed));
if (!hassplit)
{
if (maxdim < nor - 1)
MESSAGE(0,("DIMENSION LIMIT EXCEEDED\n"));
else
MESSAGE(0,("SPANS WHOLE SPACE\n"));
}
else
{
MESSAGE(0,("SUB=%6ld QUOT=%6ld\n",dim,nor-dim));
if (output)
wrmod(iseed);
}
}
else
{ if (hassplit)
{
if (smode != READALL)
MESSAGE(0,("VECTOR %ld SPLITS: ",iseed));
MESSAGE(0,("SUBSPACE %ld QUOTIENT %ld\n",dim,nor-dim));
if (output)
do_split();
exit(EXIT_OK);
}
}
}

/* ------------------------------------------------------------------
spin_up() - Spin up the seed vector (s).

returns 1 if a proper subspace was found, 0 otherwise
------------------------------------------------------------------ */

int spin_up()

{ int g;
long j, nmu, pv;
PTR xput, xget, x8;
FEL mk;

xput = m3;
dim = 0;

/* If smode == READALL, echelonize the seed space.
----------------------------------------------- */
if (smode == READALL)
{
while (dim < ssdim)
{
x8 = m3;
for (j = 1; j <= dim; ++j)
{
mk = zextract(xput,piv[(int)j]);
if (mk != F_ZERO)
zaddmulrow(xput,x8,zsub(F_ZERO,mk));
zadvance(&x8,(long)1);
}
if ((piv[(int) ++dim] = zfindpiv(xput,&mk)) == 0)
err('0');
if (mk != F_ONE)
zmulrow(xput,zdiv(F_ONE,mk));
zadvance(&xput,(long)1);
}
MESSAGE(0,("SEED HAS DIMENSION %ld\n",dim));
dim = ssdim - 1;
xput = m3;
zadvance(&xput,ssdim - 1);
}

/* Spin up
------- */
g = -1;
xget = m3;
nmu = 0;
while (1)
{
pv = zfindpiv(xput,&mk);
if (pv != 0)
{
if (mk != F_ONE)
zmulrow(xput,zdiv(F_ONE,mk));
if (dim % 100 == 0)
MESSAGE(1,("dim=%ld\n",dim));
piv[(int)++dim] = pv;
zadvance(&xput,(long)1);
if (dim > maxdim)
return 0; /* Spans whole Space */
}

/* Apply the next generator
------------------------ */
if (++g >= ngen)
{
g = 0;
zadvance(&xget,(long)1);
if (++nmu >= dim)
{
if (dim == 0) err('0');
return 1; /* Split! */
}
}
switch (gentype)
{
case MATRIX:
zmaprow(xget,gen[g],nor,xput);
break;
case PERMUTATION:
zpermrow(xget,gen[g],xput);
break;
case DONTKNOW:
err('b');
break;
}

/* Clean the image with the basis.
------------------------------- */
x8 = m3;
for (j = 1; j <= dim; ++j)
{ pv = piv[(int)j];
if ((mk = zextract(xput,pv)) != F_ZERO)
zaddmulrow(xput,x8,zsub(F_ZERO,mk));
zadvance(&x8,(long)1);
}
}
}

/* ------------------------------------------------------------------
do_split() - Calculate the action of g1 and g2 on subspace and
quotient.
------------------------------------------------------------------ */

void do_split()

{
FEL f;
int g, i;
long j, k, qdim;
PTR xa, xb;
FILE *ofile;

/* Calculate the action on the subspace
------------------------------------ */
for (g = 0; g < nmat; ++g)
{
if ((ofile = zwritehdr(subname[g],fl,dim,dim)) == NULL)
errexit(-1,subname[g]);
xa = m3;
for (j = 1; j <= dim; ++j)
{
zsetlen(fl,dim);
zmulrow(tmp2,F_ZERO);

/* Map the j-th basis vector
------------------------- */
zsetlen(fl,nor);
switch (gentype)
{
case MATRIX:
zmaprow(xa,gen[g],nor,tmp1);
break;
case PERMUTATION:
zpermrow(xa,gen[g],tmp1);
break;
case DONTKNOW:
err('b');
}
zadvance(&xa,(long)1);

/* Clean it with the basis and store the
coefficients in tmp1.
------------------------- */
for (xb = m3,k = 1; k <= dim;
zadvance(&xb,(long)1),++k)
{
f = zextract(tmp1,piv[(int)k]);
zinsert(tmp2,k,f);
if (f != F_ZERO)
zaddmulrow(tmp1,xb,zsub(F_ZERO,f));
}

/* Write one row of output
----------------------- */
zsetlen(fl,dim);
zwritevec(ofile,tmp2,(long)1);
}
fclose(ofile);
}

/* We cannot calculate the quotient for permutations
------------------------------------------------- */
if (gentype == PERMUTATION) return;

/* Find insignificant rows and extend the pivot table
-------------------------------------------------- */
qdim = nor - dim; /* Dimension of the quotient */
for (i = 1; i <= nor; ++i)
ins[i] = 1;
for (i = 1; i <= dim; ++i)
ins[(int)(piv[i])] = 0;
for (j = dim, i = 1; i <= nor; ++i)
if (ins[i])
piv[(int)++j] = i;

/* Calculate the action on the quotient
------------------------------------ */
for (g = 0; g < nmat; ++g)
{
if ((ofile = zwritehdr(quotname[g],fl,qdim,qdim)) == NULL)
errexit(-1,quotname[g]);
for (j = 1; j <= qdim; ++j)
{ zsetlen(fl,qdim);
zmulrow(tmp2,F_ZERO);

/* Take the piv[j+dim]-th row of the matrix
---------------------------------------- */
zsetlen(fl,nor);
xa = gen[g];
zadvance(&xa,piv[(int)(j+dim)]-1);
zmoverow(tmp1,xa);

/* Clean it with the subspace basis
-------------------------------- */
xb = m3;
for (k = 1; k <= dim; ++k)
{ f = zextract(tmp1,piv[(int)k]);
if (f != F_ZERO)
zaddmulrow(tmp1,xb,zsub(F_ZERO,f));
zadvance(&xb,(long)1);
}

/* Extract the insignificant columns
and write out the row.
--------------------------------- */
for (k = 1; k <= qdim; ++k)
{ f = zextract(tmp1,piv[(int)(dim+k)]);
zinsert(tmp2,k,f);
}
zsetlen(fl,qdim);
zwritevec(ofile,tmp2,(long)1);
}
fclose(ofile);
}
zsetlen(fl,nor);
}

/* ------------------------------------------------------------------
wrmod() - Write out the submodule
------------------------------------------------------------------ */

void wrmod(cnt)
long cnt;

{
char buf[200];
FILE *ofile;

if (smode == READALL)
strcpy(buf,vname);
else
{
sprintf(buf,"%s.%4.4ld",vname,cnt);
}
if ((ofile = zwritehdr(buf,fl,dim,nor)) == NULL)
errexit(-1,buf);
zwritevec(ofile,m3,dim);
fclose(ofile);
}}}