invar_space.c

#include"typedef.h"
#include"matrix.h"
#include"bravais.h"
#include"longtools.h"
/**************************************************************************\
@---------------------------------------------------------------------------
@---------------------------------------------------------------------------
@ FILE: invar_space.c
@---------------------------------------------------------------------------
@---------------------------------------------------------------------------
@
\**************************************************************************/

/**************************************************************************\
@---------------------------------------------------------------------------
@ matrix_TYP **invar_space(B, Banz, fodim, symm_opt, epsilon, anz)
@ matrix_TYP **B;
@ int Banz, fodim, symm_opt, epsilon, *anz;
@
@ calculates a basis of the lattices of integral matrices X 
@ with B[i]^{tr} * X * B[i] = X for 0<=i<Banz
@ 
@ the number of basis elements must be given to the function by the
@ integer 'fodim'
@ Then 'fodim'+1 random elements of the lattice are calculated with
@ the function 'rform' and
@ a basis of the integral matrices in the subspace generated by these
@ elements is constructed with 'long_rein_mat'.
@ This basis is the result and the number of elements is returned
@ via (int *anz)
@ The argument 'epsilon' is the only needed argument for 'rform', that
@ is called in the function.
@ the number of basis elements is returned via (int *anz).
@ 
@  'symm_opt' must be 0, 1 or -1.
@  symm_opt = 0:  the full lattice is calculated.
@  symm_opt = 1:  a basis for the symmetric matrices is calculated
@  symm_opt =-1:  a basis for the skewsymmetric matrices is calculated
@---------------------------------------------------------------------------
@
\**************************************************************************/


matrix_TYP **invar_space(B, Banz, fodim, symm_opt, epsilon, anz)
matrix_TYP **B;
int Banz, fodim, symm_opt, epsilon, *anz;
{
   int i,j,k,l, dim, ff, dd;
   matrix_TYP *pm, *pm1, **erg, *startmat, *invar;
   int r;
 
   dim = B[0]->cols;
   if(B[0]->rows != dim)
   {
     printf("error in formspace: non-square matrix in group 'B'\n");
     exit(3);
   }
   for(i=1;i<Banz;i++)
   {
      if(B[i]->rows != dim || B[i]->cols != dim)
      {
        printf("error in formspace: different dimesion of group elements\n");
        exit(3);
      }
   }
  if(symm_opt == 1)
    dd = (dim *(dim+1))/2;
  if(symm_opt == -1)
    dd = (dim *(dim-1))/2;
  if(symm_opt == 0)
    dd = dim * dim;

   if(fodim <= 0)
   { *anz = 0; return(NULL);}
  
   if(fodim == 1)
   {
     if( (erg = (matrix_TYP **)malloc(1 *sizeof(matrix_TYP *))) == NULL)
     {
       printf("malloc of 'erg' in 'invar_space' failed\n");
       exit(2);
     }
     erg[0] = init_mat(dim,dim,"");
     startmat = init_mat(dim,dim,"");
     if(symm_opt == -1 && dim > 1)
     {  startmat->array.SZ[0][1] = -1;
        startmat->array.SZ[1][0] = 1;
     }
     else
       startmat->array.SZ[0][0] = 1;
     erg[0] = rform(B, Banz, startmat, epsilon);
     free_mat(startmat);
     *anz = 1;
     return(erg);
   }

   pm = init_mat(fodim+1, dd, "");
   startmat = init_mat(dim, dim, "");
   for(i=0;i<fodim+1;i++)
   {
      if(symm_opt == 0)
      {
         for(j=0;j<dim;j++)
           for(k=0;k<dim;k++)
           {
             r = rand();
             r %= 2;
             startmat->array.SZ[j][k] = r;
           }
      }
      if(symm_opt == 1)
      {
         for(j=0;j<dim;j++)
           /* tilman: changed 14/5/97 from 
              for(k=0;k<=i;k++) to: */
           for(k=0;k<=j;k++)
           {
             r = rand();
             r %= 2;
             startmat->array.SZ[j][k] = r;
             startmat->array.SZ[k][j] = r;
           }
      }
      if(symm_opt == -1)
      {
         for(j=0;j<dim;j++)
           for(k=0;k<i;k++)
           {
             r = rand();
             r %= 2;
             startmat->array.SZ[j][k] = r;
             startmat->array.SZ[k][j] = -r;
           }
         for(j=0;j<dim;j++)
           startmat->array.SZ[j][j] = 0;
      }
      invar = rform(B, Banz, startmat, epsilon);
      if(symm_opt == 0)
      {
        l = 0;
        for(j=0;j<dim;j++)
         for(k=0;k<dim;k++)
         {  pm->array.SZ[i][l] = invar->array.SZ[j][k]; l++;}
      }
      if(symm_opt == 1)
      {
        l = 0;
        for(j=0;j<dim;j++)
         for(k=0;k<=j;k++)
         {  pm->array.SZ[i][l] = invar->array.SZ[j][k]; l++;}
      }

      /* changed 18/12/96 tilman from:
      if(symm_opt == 1)
      to: */
      if(symm_opt == (-1))
      {
        l = 0;
        for(j=0;j<dim;j++)
         for(k=0;k<j;k++)
         {  pm->array.SZ[i][l] = invar->array.SZ[j][k]; l++;}
      }
      free_mat(invar);
   }
   free_mat(startmat);

   /* output for debugging purposes
   put_mat(pm,NULL,"in invar_space, pm",2); */

   pm1 = long_rein_mat(pm);

   /* output for debugging purposes
   put_mat(pm1,NULL,"in invar_space, pm1",2); */

   free_mat(pm);

   /* changed 18/12/96 tilman from
   ff = pm->rows;
   to */
   ff = pm1->rows;

   *anz = ff;
   if( (erg = (matrix_TYP **)malloc(ff *sizeof(matrix_TYP *))) == NULL)
   {
     printf("malloc of 'erg' in 'invar_space' failed\n");
     exit(2);
   }
   for(i=0; i<ff;i++)
   {
      erg[i] = init_mat(dim,dim,"");
      l=0;
      if(symm_opt == 0)
      {
         for(j=0;j<dim;j++)
           for(k=0;k<dim;k++)
           { erg[i]->array.SZ[j][k] = pm1->array.SZ[i][l]; l++;}
      }
      if(symm_opt == 1)
      {
         for(j=0;j<dim;j++)
           for(k=0;k<=j;k++)
           {
             erg[i]->array.SZ[j][k] = pm1->array.SZ[i][l];
             erg[i]->array.SZ[k][j] = pm1->array.SZ[i][l];
             l++;
           }
      }
      if(symm_opt == -1)
      {
         for(j=0;j<dim;j++)
           for(k=0;k<j;k++)
           {
             erg[i]->array.SZ[j][k] = pm1->array.SZ[i][l];
             erg[i]->array.SZ[k][j] = pm1->array.SZ[i][l];
             l++;
           }
         for(j=0;j<dim;j++)
             erg[i]->array.SZ[j][j] = 0;
      }
      Check_mat(erg[i]);
   }
  for(i=0;i<ff;i++)
  {
    for(j=0;j<dim && erg[i]->array.SZ[j][j] == 0; j++);
    if(j<dim && erg[i]->array.SZ[j][j] < 0)
    {
       for(k=0;k<dim;k++)
        for(l=0;l<dim;l++)
          erg[i]->array.SZ[k][l] = -erg[i]->array.SZ[k][l];
    }
  }
   free_mat(pm1);
   return(erg);
}