adj.cpp

   1 // Adjunct of a matrix
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_adj(void)
   8 {
   9         push(cadr(p1));
  10         eval();
  11         adj();
  12 }
  13
  14 void
  15 adj(void)
  16 {
  17         int i, j, n;
  18
  19         save();
  20
  21         p1 = pop();
  22
  23         if (istensor(p1) && p1->u.tensor->ndim == 2 && p1->u.tensor->dim[0] == p1->u.tensor->dim[1])
  24                 ;
  25         else
  26                 stop("adj: square matrix expected");
  27
  28         n = p1->u.tensor->dim[0];
  29
  30         p2 = alloc_tensor(n * n);
  31
  32         p2->u.tensor->ndim = 2;
  33         p2->u.tensor->dim[0] = n;
  34         p2->u.tensor->dim[1] = n;
  35
  36         for (i = 0; i < n; i++)
  37                 for (j = 0; j < n; j++) {
  38                         cofactor(p1, n, i, j);
  39                         p2->u.tensor->elem[n * j + i] = pop(); /* transpose */
  40                 }
  41
  42         push(p2);
  43
  44         restore();
  45 }