CoCalc -- construct

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Matrix / construct_mat.c
⁵⁶³⁶⁸⁰ views
1
#include "typedef.h"
2
#include "matrix.h"
3
#include "tools.h"
4

5
/**************************************************************************\
6
@---------------------------------------------------------------------------
7
@---------------------------------------------------------------------------
8
@ FILE: construct_mat.c
9
@---------------------------------------------------------------------------
10
@---------------------------------------------------------------------------
11
@
12
\**************************************************************************/
13
/*
14
 |
15
 | matrices/construct.c -- Funktionen, die Matrizen erzeugen
16
 |
17
 |
18
 | exportiert die Funktionen
19
 |
20
 | init_mat
21
 | copy_mat
22
 | free_mat
23
 | Check_mat
24
 |
25
 | importiert die Funktionen
26
 | GGT, calloc2dim, malloc2dim, memcpy2dim, free2dim
27
 | importiert die Variable
28
 | act_prime
29
 |
30
*/
31

32
extern int GGT (int, int);
33
extern int act_prime;
34

35
/*{{{}}}*/
36
/*{{{  init_mat*/
37
/*{{{  Documentation*/
38
/*
39
 @
40
 @ void init_mat( cols, rows, Optionen );
41
 @ generaates a matrix with the options described below
42
 | Erzeugt eine Matrix gemaess der uebergebenen Parameter, s.u.
43
 @
44
 @ Args:
45
 @ int cols, int rows -- number of rows and columns
46
 @ char *Optionen     -- String consistin of numbers and letters 'c', 'e', 'd',
47
 @                       's', 'r', 'p'
48
 | Die Bedeutung der einzelnen Optionen:
49
 @ Explanation of the options:
50
 @
51
 @                Options
52
 @  ==========================================================================
53
 @  's'     : symmetric matrix                    |  mat->flags.Symmetric=1
54
 @  --------------------------------------------------------------------------
55
 @  'd'     : diagonal matrix, implies 's'        |  mat->flags.Diagonal=1
56
 @                                                |  mat->flags.Symmetric=1
57
 @  --------------------------------------------------------------------------
58
 @  'c', 'e': scalar matrix, implies 'd' und 's'  |  mat->flags.Scalar=1
59
 @                                                |  mat->flags.Diagonal=1
60
 @                                                |  mat->flags.Symmetric=1
61
 @  --------------------------------------------------------------------------
62
 @     numerical letter from   0 .. 9.            |  mat->flags.Scalar=1
63
 @   The same as 'c' and 'e', but the matrix      |  mat->flags.Diagonal=1
64
 @   is initialized with the number               |  mat->flags.Symmetric=1
65
 @  --------------------------------------------------------------------------
66
 @  'p'     : Matrix over  GF(p) | mat->flags.Integral= 1
67
 @                               | mat->prime= act_prime (glob. Variable)
68
 @                               | mat->array.N= NULL
69
 @  --------------------------------------------------------------------------
70
 @  'r': rational  matrix,      | mat->flags.Integral= 0
71
 @                              | mat->array.N= malloc( was auch immer );
72
 @
73
 |  falls 'r' und 'p' gleichzeitig gesetzt werden, so hat 'p' Vorrang.
74
 |  mat->flags.Integral ist immer gesetzt, es sei denn, 'r' waere
75
 |  angegeben worden oder mat->kgv != 1
76
 @  if 'r' and 'p' are used at the same time, 'p' has priority
77
 @
78
 |  mat->kgv wird bei einer Bruchmatrix als Hauptnenner aller Eintraege
79
 |  benutzt. Bei einer Matrix ueber Z ist mat->kgv= 1. Falls mat->kgv != 1,
80
 |  so ist mat->array.N = NULL.
81
 @  mat->kgv is used by a rational matrix a least common divisor
82
 @  over the denomonatop of all entries. An integral matrix has mat->kgv = 1.
83
 @  If mat->kgv != 1 then mat->array.N = NULL.
84
 @
85
 |  Nullmatrizen werden nicht besonders gekennzeichnet. Das verschwendet
86
 |  zwar Rechenzeit, ist aber Programmtechnisch weniger aufwendig.
87
 |  Insgesamt sollte bei 6x6 Matrizen die Performance nicht allzusehr
88
 |  leiden.
89
 @  Zero matrices have no extra labeling
90
 @
91
 |  Wichtig: falls im Optionen-String keine Zahlen vorkommen, so erzeugt
92
 |  init_mat() eine Nullmatrix, alle Eintraege sind also mit '0' initialisiert
93
 @  if no numerical numbers are conatained in the string, all entries 
94
 @  of the matrix are initialized with 0.
95
 @
96
*/
97
/*}}}  */
98
/*{{{  Source-Code*/
99
matrix_TYP *init_mat(rows, cols, option)
100
int cols, rows;
101
char *option;
102
{
103
int i;
104
matrix_TYP *mat;
105
int val;
106
char *temp;
107
char *tempN;
108

109
  mat = NULL;
110
  mat = (matrix_TYP *)malloc(sizeof(matrix_TYP));
111

112
  /*{{{  parse the options*/
113
  if (strpbrk(option,"pP") != NULL)
114
  {
115
    mat->prime = act_prime; /* ist auf -1 gesetzt vor dem ersten Aufruf von */
116
                            /* init_prime().                                */
117
    mat->flags.Integral = 1;
118
  }
119
  else 
120
  { 
121
    mat->prime= 0;
122
    mat->flags.Integral  = strpbrk(option, "rR") == NULL ? TRUE: FALSE;
123
  }
124
  /* changed tilman 29/07/97 to cope with negative integers from
125
  if (    (strpbrk(option, "cC") != NULL) || (strpbrk(option, "eE") != NULL)
126
       || ( (temp = strpbrk(option, "0123456789")) != NULL )
127
     ) to : */
128
  if (    (strpbrk(option, "cC") != NULL) || (strpbrk(option, "eE") != NULL)
129
       || ( (temp = strpbrk(option, "-0123456789")) != NULL )
130
     )
131
  {
132
    mat->flags.Scalar=
133
    mat->flags.Diagonal=
134
    mat->flags.Symmetric= TRUE;
135
  }
136
  else
137
  {
138
    mat->flags.Scalar= FALSE;
139
    if ( strpbrk(option, "dD") != NULL )
140
    {
141
      mat->flags.Diagonal=
142
      mat->flags.Symmetric= TRUE;
143
    } 
144
    else
145
    {
146
      mat->flags.Diagonal= FALSE;
147
      mat->flags.Symmetric = strpbrk(option, "sS") != NULL ? TRUE: FALSE;
148
    }
149
  }
150
  if ( ( tempN = strpbrk(option, "/")) != NULL ) {
151
    if ( ( tempN = strpbrk(tempN, "0123456789")) != NULL ) {
152
      mat->flags.Integral = FALSE;
153
    }
154
  }
155
  /*}}}  */
156

157
  mat->kgv= 1;
158
  mat->cols = cols;
159
  mat->rows = rows;
160
  
161
  /*{{{  alloc SZ (is always done)*/
162
  mat->array.SZ= (int **) calloc2dim( rows, cols, sizeof(int) );
163
  if(temp != NULL) 
164
  {
165
    sscanf(temp,"%d", &val);
166
    for(i = 0; i < rows; i++) mat->array.SZ[i][i] = val;
167
  }
168
  /*}}}  */
169
  if ( !mat->flags.Integral )
170
  { 
171
    if (tempN != NULL) {
172
      sscanf(tempN,"%d", &val);
173
    } else {
174
      val = 1;
175
    }
176
    mat->array.N= (int **) calloc2dim( rows, cols, sizeof(int) );
177
    memset2dim( (char **)mat->array.N,rows,cols,sizeof(int), (char *)&val);
178
  }
179
  else 
180
    mat->array.N = NULL;
181
  
182
  return(mat);
183
}
184

185
/*}}}  */
186
/*}}}  */
187
/*{{{  copy_mat*/
188
/*
189
 @-----------------------------------------------------------------
190
 @ matrix_TYP *copy_mat( matrix_TYP *old );
191
 @
192
 | kopiert die Matrix "old", Rueckgabewert ist die Kopie
193
 @ copies the matrix 'old' and returns the copy.
194
 @-----------------------------------------------------------------
195
*/
196

197
matrix_TYP *copy_mat( old )
198
matrix_TYP *old;
199
{ 
200
matrix_TYP *new= NULL;
201

202
  if ( old )
203
  { 
204
    new= (matrix_TYP *)malloc( sizeof(matrix_TYP));
205
    if ( new )
206
    { 
207
      memcpy( (char *)new, (char *)old, sizeof(matrix_TYP) );
208
      if ( new->array.SZ )
209
      { 
210
        new->array.SZ= (int **)malloc2dim( old->rows, old->cols, sizeof(int) );
211
        memcpy2dim( (char **)new->array.SZ, (char **)old->array.SZ,
212
                    old->rows, old->cols, sizeof(int) );
213
      }
214
      if ( new->array.N )
215
      { 
216
        new->array.N= (int **)malloc2dim( old->rows, old->cols, sizeof(int) );
217
        memcpy2dim( (char **)new->array.N,(char **)old->array.N,
218
                    old->rows, old->cols, sizeof(int) );
219
      }
220
    }
221
  }
222
  return new;
223
}
224

225
/*}}}  */
226
/*{{{  free_mat*/
227
/*
228
 @--------------------------------------------------------------
229
 @ void free_mat ( mat )
230
 @ matrix_TYP *mat;
231
 | -- gibt den Speicher frei, den mat benutzt
232
 @
233
 @ clear the storage allocated for mat
234
 @--------------------------------------------------------------
235
*/
236

237
void free_mat (mat)
238
matrix_TYP *mat;
239
{  
240

241
  if ( mat->array.N )
242
    free2dim( (char **)mat->array.N, mat->rows);
243
  if ( mat->array.SZ )
244
    free2dim( (char **)mat->array.SZ, mat->rows);
245
  free( (int *)mat );
246

247
}
248

249
/*}}}  */
250
/*{{{  Check_mat*/
251
/*
252
 @--------------------------------------------------------------
253
 @ void Check_mat( mat )
254
 @ matrix_TYP *mat;
255
 @
256
 | Ueberprueft die mat->flags und korrigiert diese ggf. Kuerzt
257
 | Bruchmatrizen (sowohl matrix->array.N als auch mat->kgv )
258
 | Normiert die Darstellung modulo mat->prime, so dass 0 <= Eintrag < prime
259
 @ Checks mat->flags and corrects it if necessary.
260
 @ Reduces rational matrices (mat->array.N and also mat->kgv)
261
 @ Normalizes modulo mat->prime, such that 0<= entry < prime.
262
 @--------------------------------------------------------------
263
*/
264

265
void Check_mat(mat)
266
matrix_TYP *mat;
267
{
268
int i,j;
269
int g;
270
int **Z;
271

272
  if (mat->kgv == 0) mat->kgv = 1;
273
  if ( mat->array.N )
274
  {     
275
    if ( mat->kgv != 1 ) { 
276
      fprintf(stderr,"Check_mat: Error: denominator matrix with lcd != 1 detected.\n");
277
      exit(3);
278
    }      
279
    /*
280
     *  Beim Kuerzen wird Integral richtig gesetzt.
281
     */
282
    mat->flags.Integral= TRUE;
283
    /*{{{  kuerzen*/
284
    for ( i=0; i < mat->rows; i++ ) {
285
      for ( j=0; j < mat->cols; j++ ) {
286
        if ( mat->array.N[i][j] == 0 ) {
287
          fprintf (stderr,"Check_mat: Error: divide by zero\n");
288
          exit (3);
289
        }
290
        if ( mat->array.SZ[i][j] == 0 ) {
291
          mat->array.N[i][j]= 1;        
292
        } else {
293
          g = GGT (mat->array.SZ[i][j], mat->array.N[i][j]);
294
          mat->array.SZ[i][j] /= g;
295
          mat->array.N[i][j] /= g;
296
          if ( mat->array.N[i][j] < 0) {
297
            mat->array.N[i][j] = -mat->array.N[i][j];
298
            mat->array.SZ[i][j] = -mat->array.SZ[i][j];
299
          }
300
        }
301
        if ( mat->array.N[i][j] != 1 ) {
302
          mat->flags.Integral = FALSE;
303
        }
304
      }
305
    }
306
    /*}}}  */
307
    /*{{{  mat->flags setzen*/
308
    if ( mat->cols != mat->rows ) {
309
      mat->flags.Symmetric= mat->flags.Diagonal= mat->flags.Scalar= FALSE;
310
    } else {
311
      mat->flags.Diagonal= mat->flags.Symmetric= TRUE;
312
      for ( i=0; i < mat->rows-1 && mat->flags.Symmetric;i++ ) {
313
        for ( j=i+1;j < mat->cols && mat->flags.Symmetric; j++ ) {
314
          if ( mat->array.SZ[i][j] != mat->array.SZ[j][i] || 
315
               mat->array.N [i][j] != mat->array.N [j][i]    ) {
316
            mat->flags.Symmetric= mat->flags.Diagonal= FALSE;
317
          } else if ( mat->flags.Diagonal ) {
318
            mat->flags.Diagonal= mat->array.SZ[i][j] == 0;
319
          }
320
        }
321
      }
322
      mat->flags.Scalar = mat->flags.Diagonal;
323
      for ( i=0;i < mat->rows-1 && mat->flags.Scalar; i++) {
324
        if (    mat->array.SZ[i][i] != mat->array.SZ[i+1][i+1]
325
             || mat->array.N [i][i] != mat->array.N [i+1][i+1]   ) {
326
          mat->flags.Scalar= FALSE;                 
327
        }
328
      }
329
    } 
330
    /*
331
     *   release array.N if all denominators are equal to 1
332
     */
333
    if ( mat->flags.Integral ) {
334
      free2dim( (char **)mat->array.N, mat->rows );
335
      mat->array.N = NULL;
336
    }
337
    /*}}}  */
338
  } else { /* hoechstens noch Hauptnenner */
339
    Z = mat->array.SZ;
340
    if ( mat->prime != 0 ) { /* Matrix ueber GF(p) */
341
      /*{{{  normalisieren 0 <= zahl < p*/
342
      mat->flags.Integral= 1;
343
      if ( mat->prime != -1 ) /* passiert, wenn init_mat( ... , "p" ) vor */
344
      {                       /* init_prime() aufgerufen wird             */
345
        for (i = 0; i < mat->rows; i++)
346
          for (j = 0; j < mat->cols; j++)
347
            if ( (Z[i][j] %= mat->prime) < 0)
348
              Z[i][j] += mat->prime;
349
      }   
350
      /*}}}  */
351
    } 
352
    /*{{{  flags setzen*/
353
    if ( mat->cols != mat->rows )
354
      mat->flags.Symmetric= mat->flags.Diagonal= mat->flags.Scalar= FALSE;
355
    else
356
    {
357
      mat->flags.Diagonal= mat->flags.Symmetric= TRUE;
358
      for ( i=0; i < mat->rows-1 && mat->flags.Symmetric;i++ )
359
        for ( j=i+1;j < mat->cols && mat->flags.Symmetric; j++ )
360
        {
361
          if ( mat->array.SZ[i][j] != mat->array.SZ[j][i] )
362
          { 
363
            mat->flags.Symmetric= mat->flags.Diagonal= FALSE;
364
          } 
365
          else if ( mat->flags.Diagonal )
366
            mat->flags.Diagonal= mat->array.SZ[i][j] == 0;
367
        }
368
      mat->flags.Scalar = mat->flags.Diagonal;
369
      for ( i=0;i < mat->rows-1 && mat->flags.Scalar; i++)
370
      {
371
        if ( mat->array.SZ[i][i] != mat->array.SZ[i+1][i+1] )
372
          mat->flags.Scalar= FALSE;
373
      }
374
    }
375
    /*}}}  */
376
    if ( mat->prime == 0 )
377
    {
378
      /*{{{  kgv > 0 machen*/
379
      if ( mat->kgv < 0 )
380
      {
381
        mat->kgv= - mat->kgv;
382
        Z = mat->array.SZ;
383
        if(mat->flags.Diagonal)
384
          for (i = 0; i < mat->rows; i++) Z[i][i] = -Z[i][i];
385
        else
386
          for (i = 0; i < mat->rows; i++)
387
            for (j = 0; j < mat->cols; j++)
388
              Z[i][j] = -Z[i][j];
389
      }
390
      /*}}}  */
391
      /*{{{  kgv gegen alle Matrixelemente kuerzen*/
392
      g= mat->kgv;
393
      if( mat->flags.Scalar == FALSE )
394
      {
395
        if ( mat->flags.Diagonal == FALSE )
396
        {
397
          if ( mat->flags.Symmetric == FALSE )
398
          {
399
            for (i=0; i < mat->rows && g != 1;i++ )
400
              for ( j=0;j < mat->cols && g != 1; j++ )
401
                g= GGT ( g, mat->array.SZ[i][j] );
402
          }
403
          else
404
          {
405
            for (i=0; i < mat->rows && g != 1;i++ )
406
              for ( j=0;j <=i  && g != 1; j++ )
407
                g= GGT ( g, mat->array.SZ[i][j] );
408
          }
409
        }
410
        else
411
         for ( i=0; i < mat->rows && g != 1;i++) g= GGT( g, mat->array.SZ[i][i] );
412
      }
413
      /* timan 11/05/99: changed from
414
      else
415
      to : (to handle this rare case of a matrix without entries) */
416
      else if (mat->rows > 0 && mat->cols > 0)
417
        g = GGT(g, mat->array.SZ[0][0]);
418
      mat->kgv /= g;
419
      if(g != 1)
420
      {
421
         if ( mat->flags.Diagonal )
422
           for ( i=0; i < mat->rows;i++) mat->array.SZ[i][i] /= g;
423
         else
424
           for (i=0; i < mat->rows;i++ )
425
             for ( j=0;j < mat->cols; j++ )
426
               mat->array.SZ[i][j] /= g;
427
       }
428
       mat->flags.Integral = (mat->kgv == 1);
429
    }
430
  }
431
}
432

433
/*}}}  */
434

435

436
Product

Resources

Company
