CoCalc -- inv

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

6

7
/**************************************************************************\
8
@---------------------------------------------------------------------------
9
@---------------------------------------------------------------------------
10
@ FILE: inv_mat.c
11
@---------------------------------------------------------------------------
12
@---------------------------------------------------------------------------
13
@
14
\**************************************************************************/
15
/*{{{}}}*/
16
/*{{{  rmat_inv       */
17
/*
18
*  result = rmat_inv( mat );
19
*
20
*  berechnet Inverse zu eine Nennermatrix. Das "r" steht fuer rational
21
*
22
*  matrix_TYP *result: Inverse zu mat, wird erzeugt
23
*  matrix_TYP *mat:    wird nicht veraendert
24
 */
25
static matrix_TYP *rmat_inv ( mat )
26
matrix_TYP *mat;                           
27
{  
28
rational **M, **II, *v, f ;
29
int n, i, j, k ,l , flag;
30
int *nuin_M, *nuin_I;
31
int kgv, h;
32
matrix_TYP *inv;
33

34
  flag  =
35
  kgv   =
36
  h     = 0;
37
  f = Zero;
38
  n = mat->cols;
39
  nuin_M = (int *)calloc((n+1),sizeof(int));
40
  nuin_I = (int *)calloc((n+1),sizeof(int));
41
  
42
  /*
43
   * create rational version II of identity matrix        
44
   */
45
  /*{{{  */
46
  II = ( rational ** )malloc2dim( n, n, sizeof(rational) );
47
  memset2dim( (char ** )II, n, n, sizeof(rational), (char *)&Zero );
48
  for(i = 0; i < n; i++) {
49
    II[i][i] = One;
50
  }
51
  /*}}}  */
52
  /*
53
   * create rational version M of matrix mat
54
   */
55
  /*{{{  */
56
  M  = (rational **)malloc2dim( n, n, sizeof(rational) );
57
  for ( i  = 0; i  < n; i ++) {
58
    for ( j  = 0; j  < n; j ++) {
59
      M[i][j].z = mat->array.SZ[i][j];
60
      M[i][j].n = mat->array.N[i][j];
61
    }
62
  }
63
  /*}}}  */
64
  /*
65
   * Invert M by paralell transformation of M and II s.t M is unit 
66
   */                     
67

68
  /*
69
   *  first create triangular matrix
70
   */
71
  /*{{{  */
72
  for (i = 0; i < n; i++) {
73
    if ( M[i][i].z == 0 ) {
74
      /*
75
       * Swap a non-zero element into [i][i]-position        
76
       */
77
      /*{{{  */
78
      for (j = i+1; j < n && M[j][i].z == 0; j++ );
79
      if ( j == n ) {
80
        fprintf (stderr, "matrix is singular. Can't invert\n");
81
        exit (3);
82
      }
83
      v   =  M[i];  M[i] =  M[j];  M[j] = v;
84
      v   = II[i]; II[i] = II[j]; II[j] = v;
85
      /*}}}  */
86
    }
87
    /*
88
     * Find the non-zero entries in i-th row in M and II
89
     * 17.02.92, stored in var "f"
90
     */
91
    /*{{{  */
92
    nuin_M[0] =
93
    nuin_I[0] = 1;
94
    for(k = 0; k < n; k++) {
95
      if ( M[i][k].z != 0) nuin_M[nuin_M[0]++] = k;
96
      if (II[i][k].z != 0) nuin_I[nuin_I[0]++] = k;
97
    }
98
    f.z= M[i][i].z;
99
    f.n= M[i][i].n;
100
    /*}}}  */
101
    if ( (f.z != 1) || (f.n != 1) ) {
102
      /*
103
       * Normalize i-th row s.t. [i][i]-element is 1
104
       */
105
      /*{{{  */
106
      for (k = nuin_M[1], l = 1; l < nuin_M[0]; k=nuin_M[++l]) {
107
        M[i][k].z *= f.n;
108
        M[i][k].n *= f.z;
109
        Normal (&M[i][k]);
110
      }
111
      for (k = nuin_I[1], l = 1; l < nuin_I[0]; k=nuin_I[++l]) {
112
        II[i][k].z *= f.n;
113
        II[i][k].n *= f.z;
114
        Normal (&II[i][k]);
115
      }                              
116
      /*}}}  */
117
    }
118
    /*
119
     * Clear i-th column downwards
120
     */
121
    /*{{{  */
122
    for ( j = i+1; j < n; j++) {
123
      if ( M[j][i].z != 0 ) {
124
        f.z = M[j][i].z;
125
        f.n = M[j][i].n;
126
        for (k = nuin_M[1], l = 1; l < nuin_M[0]; k = nuin_M[++l]) {
127
          if(M[j][k].z != 0) {
128
            h= f.n * M[i][k].n;
129
            M[j][k].z *= h;
130
            M[j][k].z -= M[j][k].n * M[i][k].z * f.z;
131
            M[j][k].n *= h;
132
          } else {
133
            M[j][k].z= - M[i][k].z * f.z;
134
            M[j][k].n = M[i][k].n * f.n;
135
          }
136
          Normal (&M[j][k]);
137
        }
138
        for (k = nuin_I[1], l = 1; l < nuin_I[0]; k = nuin_I[++l]) {
139
          if(II[j][k].z != 0) {
140
            h= f.n * II[i][k].n;
141
            II[j][k].z *= h;
142
            II[j][k].z -= II[j][k].n * II[i][k].z * f.z;
143
            II[j][k].n *= h;
144
          } else {
145
            II[j][k].z = - II[i][k].z*f.z;
146
            II[j][k].n = II[i][k].n * f.n;
147
          }
148
          Normal (&II[j][k]);
149
        }
150
      }
151
    }
152
    /*}}}  */
153
  } /* M is now triangular */
154
  /*}}}  */
155
  /*
156
   * Clear columns upwards  
157
   */
158
  /*{{{  */
159
  for (i = n-1; i >= 0; i --) {
160
    nuin_I[0] = 1;
161
    for ( k = 0; k < n; k++) {
162
      if (II[i][k].z != 0) nuin_I[nuin_I[0]++] = k;
163
    }
164
    for ( j = i-1; j >= 0; j--) {
165
      if ( M[j][i].z != 0 ) {
166
        for (k = nuin_I[1], l = 1; l < nuin_I[0]; k = nuin_I[++l]) {
167
          h= M[j][i].n* II[i][k].n;
168
          II[j][k].z *= h;
169
          II[j][k].z -= M[j][i].z * II[j][k].n * II[i][k].z;
170
          II[j][k].n *= h;
171
          Normal (&II[j][k]);
172
        }                        
173
      }
174
    }
175
  }
176
  /*}}}  */
177
  /*
178
   * copy M to result matrix "inv"
179
   *
180
   */
181
  inv = init_mat(n,n,"r");
182
  for ( i  = 0; i < n  ; i++) {
183
    for ( j = 0; j < n; j++) {
184
      inv->array.SZ[i][j] = II[i][j].z;
185
      inv->array.N [i][j] = II[i][j].n;
186
    }
187
  }
188
  free2dim( (char **)M, n );
189
  free2dim( (char **)II, n );
190
  free(nuin_M);
191
  free(nuin_I);
192
  
193
  Check_mat( inv );
194

195
  return ( inv );
196
}
197

198
/*}}}  */
199
#ifdef __OLD__
200
/*{{{  imat_inv */
201
static matrix_TYP *imat_inv (mat)
202
matrix_TYP *mat;
203
{  
204
rational **M, **II, *v, f ;
205
int n, i, j, k ,l , flag;
206
int *nuin_M, *nuin_I;
207
int kgv, h, waste ;
208
matrix_TYP *inv;
209

210
  flag  =
211
  kgv   =
212
  h     =
213
  waste = 0;
214
  f = Zero;
215
  n = mat->cols;
216
  nuin_M = (int *)calloc((n+1),sizeof(int),"imat_inv:nuin_M");
217
  nuin_I = (int *)calloc((n+1),sizeof(int),"imat_inv:nuin_I");
218
  
219
  /*
220
   * create rational version II of identity matrix        
221
   */
222
  /*{{{  */
223
  II    = (rational **)malloc(n*sizeof(rational *),"imat_inv:II");
224
  II[0] = (rational *)malloc(n*sizeof(rational),"imat_inv:II[0]");
225
  for (i = 0; i < n; i++) {
226
    II[0][i] = Zero;
227
  }
228
  for(i = 1; i < n; i++) {
229
    II[i] = (rational *)malloc(n*sizeof(rational),"imat_inv:II[i]");
230
    memcpy(II[i],II[0],n*sizeof(rational));
231
    II[i][i] = One;
232
  }
233
  II[0][0] = One;
234

235
  /*}}}  */
236
  /*
237
   * create rational version M of matrix mat
238
   */
239
  /*{{{  */
240
  M  = (rational **)malloc(n*sizeof(rational *),"imat_inv:M");
241
  if(mat->flags.Integral) {
242
    for ( i  = 0; i  < n; i ++) {
243
      M[i]  = (rational *)malloc(n*sizeof(rational),"imat_inv:M[i]");
244
      for ( j  = 0; j  < n; j ++) {
245
        M[i][j].z = mat->array.SZ[i][j];
246
        M[i][j].n = 1;
247
      }
248
    }
249
  } else {
250
    for ( i  = 0; i  < n; i ++) {
251
      M[i]  = (rational *)malloc(n*sizeof(rational),"imat_inv:M[i]");
252
      for ( j  = 0; j  < n; j ++) {
253
        M[i][j].z = mat->array.SZ[i][j];
254
        if ( M[i][j].z == 0 ) {
255
          M[i][j].n = 1;
256
        } else {
257
          M[i][j].n = mat->kgv;
258
          if ( abs(M[i][j].z) > 1) {
259
            Normal(&M[i][j]);
260
          }
261
        }
262
      }
263
    }
264
  }
265
  /*}}}  */
266
  /*
267
   * Invert M by paralell transformation of M and II s.t M is unit 
268
   */                     
269

270
  /*
271
   *  first create triangular matrix
272
   */
273
  /*{{{  */
274
  for (i = 0; i < n; i++) {
275
    if ( M[i][i].z == 0 ) {
276
      /*
277
       * Swap a non-zero element into [i][i]-position        
278
       */
279
      /*{{{  */
280
      for (j = i+1; j < n && M[j][i].z != 0; j++ );
281
      if ( j == n ) {
282
        fprintf (stderr, "matrix is singular. Can't invert\n");
283
        exit (3);
284
      }
285
      v   =  M[i];  M[i] =  M[j];  M[j] = v;
286
      v   = II[i]; II[i] = II[j]; II[j] = v;
287
      /*}}}  */
288
    }
289
    /*
290
     * Find the non-zero entries in i-th row in M and II
291
     * 17.02.92, stored in var "f"
292
     */
293
    /*{{{  */
294
    nuin_M[0] =
295
    nuin_I[0] = 1;
296
    for(k = 0; k < n; k++) {
297
      if ( M[i][k].z != 0) nuin_M[nuin_M[0]++] = k;
298
      if (II[i][k].z != 0) nuin_I[nuin_I[0]++] = k;
299
    }
300
    f.z= M[i][i].z;
301
    f.n= M[i][i].n;
302
    /*}}}  */
303
    if ( (f.z != 1) || (f.n != 1) ) {
304
      /*
305
       * Normalize i-th row s.t. [i][i]-element is 1
306
       */
307
      /*{{{  */
308
      for (k = nuin_M[1], l = 1; l < nuin_M[0]; k=nuin_M[++l]) {
309
        M[i][k].z *= f.n;
310
        M[i][k].n *= f.z;
311
        Normal (&M[i][k]);
312
      }
313
      for (k = nuin_I[1], l = 1; l < nuin_I[0]; k=nuin_I[++l]) {
314
        II[i][k].z *= f.n;
315
        II[i][k].n *= f.z;
316
        Normal (&II[i][k]);
317
      }                              
318
      /*}}}  */
319
    }
320
    /*
321
     * Clear i-th column downwards
322
     */
323
    /*{{{  */
324
    for ( j = i+1; j < n; j++) {
325
      if ( M[j][i].z != 0 ) {
326
        f.z = M[j][i].z;
327
        f.n = M[j][i].n;
328
        for (k = nuin_M[1], l = 1; l < nuin_M[0]; k = nuin_M[++l]) {
329
          if(M[j][k].z != 0) {
330
            h= f.n * M[i][k].n;
331
            M[j][k].z *= h;
332
            M[j][k].z -= M[j][k].n * M[i][k].z * f.z;
333
            M[j][k].n *= h;
334
          } else {
335
            M[j][k].z= - M[i][k].z * f.z;
336
            M[j][k].n = M[i][k].n * f.n;
337
          }
338
          Normal (&M[j][k]);
339
        }
340
        for (k = nuin_I[1], l = 1; l < nuin_I[0]; k = nuin_I[++l]) {
341
          if(II[j][k].z != 0) {
342
            h= f.n * II[i][k].n;
343
            II[j][k].z *= h;
344
            II[j][k].z -= II[j][k].n * II[i][k].z * f.z;
345
            II[j][k].n *= h;
346
          } else {
347
            II[j][k].z = - II[i][k].z*f.z;
348
            II[j][k].n = II[i][k].n * f.n;
349
          }
350
          Normal (&II[j][k]);
351
        }
352
      }
353
    }
354
    /*}}}  */
355
  } /* end of inversion loop */
356
  /*}}}  */
357
  /*
358
   * Clear columns upwards  and calculate lcm of denominators    
359
   */
360
  kgv= II[n-1][0].n;
361
  for (i = n-1; i >= 0; i --) {
362
    nuin_I[0] = 1;
363
    for ( k = 0; k < n; k++) {
364
      kgv = kgv / GGT(kgv, II[i][k].n) * II[i][k].n;
365
      if (II[i][k].z != 0) nuin_I[nuin_I[0]++] = k;
366
    }
367
    for ( j = i-1; j >= 0; j--) {
368
      if ( M[j][i].z != 0 ) {
369
        for (k = nuin_I[1], l = 1; l < nuin_I[0]; k = nuin_I[++l]) {
370
          h= M[j][i].n* II[i][k].n;
371
          II[j][k].z *= h;
372
          II[j][k].z -= M[j][i].z * II[j][k].n * II[i][k].z;
373
          II[j][k].n *= h;
374
          Normal (&II[j][k]);
375
        }                        
376
      }
377
    }
378
  }
379
  
380
  inv = init_mat(n,n,"ik");
381
    /*------------------------------------------------------------*\
382
    | Normalize all entries to the form II[i][j].z / kgv      |
383
    | and copy into inv->array.SZ                  |
384
    \*------------------------------------------------------------*/
385
  for ( i  = 0; i < n  ; i++) {
386
    for ( j = 0; j < n; j++) {
387
      II[i][j].z = II[i][j].z * kgv /II[i][j].n;
388
      inv->array.SZ[i][j] = II[i][j].z;
389
    }
390
    free(M[i]);
391
    free(II[i]);
392
  }
393
  free(M);
394
  free(II);
395
  free(nuin_M);
396
  free(nuin_I);
397
  
398
  inv->kgv     =  kgv;
399
  inv->flags.Integral = (kgv == 1);
400
  inv->prime = mat->prime;
401
  inv->Type  = inv->Type & ~(I_Typ | P_Typ);
402
  Check_mat( inv );
403
  return ( inv );
404
}
405

406
/*}}}  */
407
#else
408
/*{{{  imat_inv*/
409
static matrix_TYP *imat_inv (mat)
410
matrix_TYP *mat;
411
{  
412
matrix_TYP *inv, *help;
413

414
  /* don't relie on kgv2rat() and rat2kgv() being the inverse of each other
415
   * there might be integer overflows and we don't want to demolish the 
416
   * original matrix
417
   */
418
  help = copy_mat( mat );
419
  /* 
420
   * Allocate a real rational denominator matrix
421
   */
422
  kgv2rat( help );
423
  /*
424
   * compute the invers of it
425
   */
426
  inv = rmat_inv( help );    
427
  /*
428
   *  free workspace
429
   */
430
  free_mat( help );
431
  /*
432
   *  transform the result to lcm-representation. This might be stupid. 
433
   *  There could be integer overflow and in that case we'd better stay to the
434
   *  denominator matrix
435
   */
436
  rat2kgv( inv );       
437
  /*
438
   *  eigentlich sollte ein Check_mat() noetig sein. rat2kgv() berechnet 
439
   *  schon eine maximal gekuerzte Darstellung
440
   */                                         
441
#if 0
442
  Check_mat( inv );
443
#endif
444
  return ( inv );
445
}
446

447
#endif
448

449
#ifdef __OLD__
450

451
/* replaced the above function to use mp integres., tilman 02/05/97 */
452
static matrix_TYP *imat_inv (mat)
453
matrix_TYP *mat;
454
{
455
   matrix_TYP *ID,
456
             **X,
457
              *res;
458

459
   ID = init_mat(mat->cols,mat->cols,"1");
460

461
   X = long_solve_mat(mat,ID);
462

463
   free_mat(ID);
464

465
   if (X == NULL){
466
      fprintf(stderr,"matrix is not invertible\n");
467
      exit(3);
468
   }
469

470
   res = X[0];
471

472
   free(X);
473

474
   return res;
475
}
476

477
/*}}}  */
478
#endif
479
/*{{{  pmat_inv, unchanged*/
480

481

482
/**************************************************************************\
483
@---------------------------------------------------------------------------
484
@ matrix_TYP *pmat_inv (mat)
485
@ matrix_TYP *mat;
486
@
487
@ calculates the inverse of mat modulo mat->prime
488
@---------------------------------------------------------------------------
489
@
490
\**************************************************************************/
491
matrix_TYP *pmat_inv (mat)
492
matrix_TYP *mat;
493
{
494
int n, f, *v;
495
int i, j, k, **M, **II;
496
matrix_TYP *inv;
497

498
  n = mat->cols;
499
  init_prime(mat->prime);
500

501
  M = (int **)malloc2dim( n, n, sizeof(int) );
502
  memcpy2dim( (char **)M, (char **)mat->array.SZ, n, n, sizeof(int) );
503
  II = (int **)calloc2dim( n, n, sizeof(int) );
504
  for ( i= 0; i< n; i ++) {
505
    II[i][i] = 1;
506
  }
507
  /*
508
   * Invert M by paralell transformation of M and II s.t M is unit 
509
   */
510
  for ( i= 0; i< n; i ++ ) {
511
    if ( M[j=i][j] == 0 ) {
512
   /*
513
    * Swap a non-zero element into [i][i]-position 
514
    */
515
      for ( j= i+1; j< n && M[j][i] == 0; j++ );
516
      if ( j == n ) {
517
        fprintf (stderr, "matrix is singular. Can't invert\n");
518
        exit (3);
519
      }
520
      v= M[i];   M[i] =  M[j];  M[j] = v;
521
      v= II[i]; II[i] = II[j]; II[j] = v;
522
    }
523
    f = M[i][i];
524
    if (f != 1) {
525
    /*
526
     * Normalize i-th row s.t. [i][i]-element is 1
527
     */
528
      for ( k= 0; k< n; k ++ ) {
529
        if ( k >= i ) {
530
          M[i][k] = P(M[i][k], - f);
531
        }
532
        II[i][k] = P(II[i][k], - f);
533
      }          
534
    }
535
    /*
536
     * Clear i-th column                     
537
     */
538
    for ( j= 0; j< n; j ++ ) {
539
      if ((j != i) && (f = M[j][i]) ) {
540
        for ( k= 0; k< n; k ++ ) {
541
          if (k > i) {
542
            M[j][k] = S(M[j][k], - P(f,M[i][k]));
543
          }
544
          II[j][k] = S(II[j][k], - P(f,II[i][k]));
545
        }
546
      }
547
    }
548
  }
549
  inv = init_mat(n,n,"p");
550
  free2dim( (char **)M, n );
551
  free2dim( (char **)inv->array.SZ, n  );
552
  inv->array.SZ = II;
553
  
554
  inv->prime = mat->prime;
555

556
  inv->flags.Symmetric = mat->flags.Symmetric;
557
  inv->flags.Diagonal  = mat->flags.Diagonal ;
558
  inv->flags.Scalar    = mat->flags.Scalar ;
559
  
560
  return (inv);
561
}
562

563
/*}}}  */
564
/*{{{  mat_inv*/
565
/**************************************************************************\
566
@---------------------------------------------------------------------------
567
@ matrix_TYP *mat_inv (mat)
568
@ matrix_TYP *mat;
569
@
570
@ calculates the inverse of mat
571
@---------------------------------------------------------------------------
572
@
573
\**************************************************************************/
574
matrix_TYP *mat_inv(mat)
575
matrix_TYP *mat;
576
{
577
matrix_TYP *inv;
578

579
  if ( mat->rows != mat->cols ) {
580
    fprintf (stderr, "Can't invert non-square matrix\n");
581
    exit (3);
582
  }
583
  if ( mat->prime > 0 ) {
584
    inv = pmat_inv( mat );
585
  } else if ( mat->array.N != NULL ) {
586
    inv = rmat_inv( mat );
587
  } else {
588
    inv = imat_inv( mat );
589
  }
590
  return(inv);
591
}
592
/*}}}  */
593

594
Product

Resources

Company
