source/kernel/matrix.c

   1 /* matrix.c
   2      COPYRIGHT
   3           Both this software and its documentation are
   4
   5               Copyright 1993 by IRISA /Universite de Rennes I -
   6               France, Copyright 1995,1996 by BYU, Provo, Utah
   7                          all rights reserved.
   8
   9 */
  10
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <ctype.h>
  15 #include <polylib/polylib.h>
  16
  17 #ifdef mac_os
  18   #define abs __abs
  19 #endif
  20
  21 /*
  22  * Allocate space for matrix dimensioned by 'NbRows X NbColumns'.
  23  */
  24 Matrix *Matrix_Alloc(unsigned NbRows,unsigned NbColumns) {
  25
  26   Matrix *Mat;
  27   Value *p, **q;
  28   int i;
  29
  30   Mat=(Matrix *)malloc(sizeof(Matrix));
  31   if(!Mat) {
  32     errormsg1("Matrix_Alloc", "outofmem", "out of memory space");
  33     return 0;
  34   }
  35   Mat->NbRows=NbRows;
  36   Mat->NbColumns=NbColumns;
  37   if (NbRows==0 || NbColumns==0) {
  38       Mat->p = (Value **)0;
  39       Mat->p_Init= (Value *)0;
  40       Mat->p_Init_size = 0;
  41   } else {
  42       q = (Value **)malloc(NbRows * sizeof(*q));
  43       if(!q) {
  44         free(Mat);
  45         errormsg1("Matrix_Alloc", "outofmem", "out of memory space");
  46         return 0;
  47       }
  48       p = value_alloc(NbRows * NbColumns, &Mat->p_Init_size);
  49       if(!p) {
  50         free(q);
  51         free(Mat);
  52         errormsg1("Matrix_Alloc", "outofmem", "out of memory space");
  53         return 0;
  54       }
  55       Mat->p = q;
  56       Mat->p_Init = p;
  57       for (i=0;i<NbRows;i++) {
  58         *q++ = p;
  59         p += NbColumns;
  60       }
  61   }
  62   p = NULL;
  63   q = NULL;
  64
  65   return Mat;
  66 } /* Matrix_Alloc */
  67
  68 /*
  69  * Free the memory space occupied by Matrix 'Mat'
  70  */
  71 void Matrix_Free(Matrix *Mat)
  72 {
  73   if (Mat->p_Init)
  74     value_free(Mat->p_Init, Mat->p_Init_size);
  75
  76   if (Mat->p)
  77     free(Mat->p);
  78   free(Mat);
  79
  80 } /* Matrix_Free */
  81
  82 void Matrix_Extend(Matrix *Mat, unsigned NbRows)
  83 {
  84   Value *p, **q;
  85   int i;
  86
  87   q = (Value **)realloc(Mat->p, NbRows * sizeof(*q));
  88   if(!q) {
  89     errormsg1("Matrix_Extend", "outofmem", "out of memory space");
  90     return;
  91   }
  92   Mat->p = q;
  93   if (Mat->p_Init_size < NbRows * Mat->NbColumns) {
  94     p = (Value *)realloc(Mat->p_Init, NbRows * Mat->NbColumns * sizeof(Value));
  95     if(!p) {
  96       errormsg1("Matrix_Extend", "outofmem", "out of memory space");
  97       return;
  98     }
  99     Mat->p_Init = p;
 100     Vector_Set(Mat->p_Init + Mat->NbRows*Mat->NbColumns, 0,
 101                Mat->p_Init_size - Mat->NbRows*Mat->NbColumns);
 102     for (i = Mat->p_Init_size; i < Mat->NbColumns*NbRows; ++i)
 103         value_init(Mat->p_Init[i]);
 104     Mat->p_Init_size = Mat->NbColumns*NbRows;
 105   } else
 106     Vector_Set(Mat->p_Init + Mat->NbRows*Mat->NbColumns, 0,
 107                (NbRows - Mat->NbRows) * Mat->NbColumns);
 108   for (i=0;i<NbRows;i++) {
 109     Mat->p[i] = Mat->p_Init + (i * Mat->NbColumns);
 110   }
 111   Mat->NbRows = NbRows;
 112 }
 113
 114 /*
 115  * Print the contents of the Matrix 'Mat'
 116  */
 117 void Matrix_Print(FILE *Dst, const char *Format, Matrix *Mat)
 118 {
 119   Value *p;
 120   int i, j;
 121   unsigned NbRows, NbColumns;
 122
 123   fprintf(Dst,"%d %d\n", NbRows=Mat->NbRows, NbColumns=Mat->NbColumns);
 124   if (NbColumns ==0) {
 125     fprintf(Dst, "\n");
 126     return;
 127   }
 128   for (i=0;i<NbRows;i++) {
 129     p=*(Mat->p+i);
 130     for (j=0;j<NbColumns;j++) {
 131       if (!Format) {
 132         value_print(Dst," "P_VALUE_FMT" ",*p++);
 133       }
 134       else {
 135         value_print(Dst,Format,*p++);
 136       }
 137     }
 138     fprintf(Dst, "\n");
 139   }
 140 } /* Matrix_Print */
 141
 142 /*
 143  * Read the contents of the Matrix 'Mat' from standard input
 144  */
 145 Matrix *Matrix_Read_Input(Matrix *Mat) {
 146   return Matrix_Read_InputFile(Mat, stdin);
 147 } /* Matrix_Read_Input */
 148
 149 /*
 150  * Read the contents of the Matrix 'Mat' from file open in fp
 151  */
 152 Matrix *Matrix_Read_InputFile(Matrix *Mat, FILE *fp) {
 153
 154   Value *p;
 155   int i,j,n;
 156   char *c, s[1024];
 157
 158   p = Mat->p_Init;
 159   for (i=0;i<Mat->NbRows;i++) {
 160     do
 161     {
 162       c = fgets(s, 1024, fp);
 163       if( !c )
 164         break;
 165       /* jump to the first non space char */
 166       while(isspace(*c) && *c!='\n' && *c)
 167         ++c;
 168     }
 169     while(*c =='#' || *c== '\n');       /* continue if the line is empty or is a comment */
 170     if (!c) {
 171       errormsg1( "Matrix_Read", "baddim", "not enough rows" );
 172       return(NULL);
 173     }
 174
 175     if( c-s >= 1023 )
 176     {
 177       /* skip to EOL and ignore the rest of the line if longer than 1024 */
 178       errormsg1( "Matrix_Read", "warning", "line longer than 1024 chars (ignored remaining chars till EOL)" );
 179       do
 180         n = getchar();
 181       while( n!=EOF && n!='\n' );
 182     }
 183
 184
 185     for (j=0;j<Mat->NbColumns;j++)
 186     {
 187       char *z;
 188       if(*c=='\n' || *c=='#' || *c=='\0') {
 189         errormsg1("Matrix_Read", "baddim", "not enough columns");
 190         return(NULL);
 191       }
 192       /* go the the next space or end */
 193       for( z=c ; *z ; z++ )
 194       {
 195         if( *z=='\n' || *z=='#' || isspace(*z) )
 196           break;
 197       }
 198       if( *z )
 199         *z = '\0';
 200       else
 201         z--; /* hit EOL :/ go back one char */
 202       value_read(*(p++),c);
 203       /* go point to the next non space char */
 204       c = z+1;
 205       while( isspace(*c) )
 206         c++;
 207     }
 208   }
 209   return( Mat );
 210 } /* Matrix_Read_InputFile */
 211
 212 /*
 213  * Read the contents of the matrix 'Mat' from standard input.
 214  * a '#' is a comment (till EOL)
 215  */
 216 Matrix *Matrix_Read(void) {
 217   return Matrix_ReadFile(stdin);
 218 } /* Matrix_Read */
 219
 220 /*
 221  * Read the contents of the matrix 'Mat' from file open in fp.
 222  * a '#' is a comment (till EOL)
 223  */
 224 Matrix *Matrix_ReadFile(FILE *fp){
 225   Matrix *Mat;
 226   unsigned NbRows, NbColumns;
 227   char s[1024];
 228
 229   if (fgets(s, 1024, fp) == NULL)
 230     return NULL;
 231   while ((*s=='#' || *s=='\n') ||
 232          (sscanf(s, "%d %d", &NbRows, &NbColumns)<2)) {
 233     if (fgets(s, 1024, fp) == NULL)
 234       return NULL;
 235   }
 236   Mat = Matrix_Alloc(NbRows,NbColumns);
 237   if(!Mat) {
 238     errormsg1("Matrix_Read", "outofmem", "out of memory space");
 239     return(NULL);
 240   }
 241   if( !Matrix_Read_Input(Mat) )
 242   {
 243     Matrix_Free(Mat);
 244     Mat=NULL;
 245   }
 246   return Mat;
 247 } /* Matrix_ReadFile */
 248
 249 /*
 250  * Basic hermite engine
 251  */
 252 static int hermite(Matrix *H,Matrix *U,Matrix *Q) {
 253
 254   int nc, nr, i, j, k, rank, reduced, pivotrow;
 255   Value pivot,x,aux;
 256   Value *temp1, *temp2;
 257
 258   /*                     T                     -1   T */
 259   /* Computes form: A = Q H  and U A = H  and U  = Q  */
 260
 261   if (!H) {
 262     errormsg1("Domlib", "nullH", "hermite: ? Null H");
 263     return -1;
 264   }
 265   nc = H->NbColumns;
 266   nr = H->NbRows;
 267   temp1 = (Value *) malloc(nc * sizeof(Value));
 268   temp2 = (Value *) malloc(nr * sizeof(Value));
 269   if (!temp1 ||!temp2) {
 270     errormsg1("Domlib", "outofmem", "out of memory space");
 271     return -1;
 272   }
 273
 274   /* Initialize all the 'Value' variables */
 275   value_init(pivot); value_init(x);
 276   value_init(aux);
 277   for(i=0;i<nc;i++)
 278     value_init(temp1[i]);
 279   for(i=0;i<nr;i++)
 280     value_init(temp2[i]);
 281
 282 #ifdef DEBUG
 283   fprintf(stderr,"Start  -----------\n");
 284   Matrix_Print(stderr,0,H);
 285 #endif
 286   for (k=0, rank=0; k<nc && rank<nr; k=k+1) {
 287     reduced = 1;        /* go through loop the first time */
 288 #ifdef DEBUG
 289     fprintf(stderr, "Working on col %d.  Rank=%d ----------\n", k+1, rank+1);
 290 #endif
 291     while (reduced) {
 292       reduced=0;
 293
 294       /* 1. find pivot row */
 295       value_absolute(pivot,H->p[rank][k]);
 296
 297       /* the kth-diagonal element */
 298       pivotrow = rank;
 299
 300       /* find the row i>rank with smallest nonzero element in col k */
 301       for (i=rank+1; i<nr; i++) {
 302         value_absolute(x,H->p[i][k]);
 303         if (value_notzero_p(x) &&
 304             (value_lt(x,pivot) || value_zero_p(pivot))) {
 305           value_assign(pivot,x);
 306           pivotrow = i;
 307         }
 308       }
 309
 310       /* 2. Bring pivot to diagonal (exchange rows pivotrow and rank) */
 311       if (pivotrow != rank) {
 312         Vector_Exchange(H->p[pivotrow],H->p[rank],nc);
 313         if (U)
 314           Vector_Exchange(U->p[pivotrow],U->p[rank],nr);
 315         if (Q)
 316           Vector_Exchange(Q->p[pivotrow],Q->p[rank],nr);
 317
 318 #ifdef DEBUG
 319         fprintf(stderr,"Exchange rows %d and %d  -----------\n", rank+1, pivotrow+1);
 320         Matrix_Print(stderr,0,H);
 321 #endif
 322       }
 323       value_assign(pivot,H->p[rank][k]);        /* actual ( no abs() ) pivot */
 324
 325       /* 3. Invert the row 'rank' if pivot is negative */
 326       if (value_neg_p(pivot)) {
 327         value_oppose(pivot,pivot); /* pivot = -pivot */
 328         for (j=0; j<nc; j++)
 329           value_oppose(H->p[rank][j],H->p[rank][j]);
 330
 331         /* H->p[rank][j] = -(H->p[rank][j]); */
 332         if (U)
 333           for (j=0; j<nr; j++)
 334             value_oppose(U->p[rank][j],U->p[rank][j]);
 335
 336         /* U->p[rank][j] = -(U->p[rank][j]); */
 337         if (Q)
 338           for (j=0; j<nr; j++)
 339             value_oppose(Q->p[rank][j],Q->p[rank][j]);
 340
 341         /* Q->p[rank][j] = -(Q->p[rank][j]); */
 342 #ifdef DEBUG
 343         fprintf(stderr,"Negate row %d  -----------\n", rank+1);
 344         Matrix_Print(stderr,0,H);
 345 #endif
 346
 347       }
 348       if (value_notzero_p(pivot)) {
 349
 350         /* 4. Reduce the column modulo the pivot */
 351         /*    This eventually zeros out everything below the */
 352         /*    diagonal and produces an upper triangular matrix */
 353
 354         for (i=rank+1;i<nr;i++) {
 355           value_assign(x,H->p[i][k]);
 356           if (value_notzero_p(x)) {
 357             value_modulus(aux,x,pivot);
 358
 359             /* floor[integer division] (corrected for neg x) */
 360             if (value_neg_p(x) && value_notzero_p(aux)) {
 361
 362               /* x=(x/pivot)-1; */
 363               value_division(x,x,pivot);
 364               value_decrement(x,x);
 365             }
 366             else
 367               value_division(x,x,pivot);
 368             for (j=0; j<nc; j++) {
 369               value_multiply(aux,x,H->p[rank][j]);
 370               value_subtract(H->p[i][j],H->p[i][j],aux);
 371             }
 372
 373             /* U->p[i][j] -= (x * U->p[rank][j]); */
 374             if (U)
 375               for (j=0; j<nr; j++) {
 376                 value_multiply(aux,x,U->p[rank][j]);
 377                 value_subtract(U->p[i][j],U->p[i][j],aux);
 378               }
 379
 380             /* Q->p[rank][j] += (x * Q->p[i][j]); */
 381             if (Q)
 382               for(j=0;j<nr;j++) {
 383                 value_addmul(Q->p[rank][j], x, Q->p[i][j]);
 384               }
 385             reduced = 1;
 386
 387 #ifdef DEBUG
 388             fprintf(stderr,
 389                     "row %d = row %d - %d row %d -----------\n", i+1, i+1, x, rank+1);
 390             Matrix_Print(stderr,0,H);
 391 #endif
 392
 393           } /* if (x) */
 394         } /* for (i) */
 395       } /* if (pivot != 0) */
 396     } /* while (reduced) */
 397
 398     /* Last finish up this column */
 399     /* 5. Make pivot column positive (above pivot row) */
 400     /*    x should be zero for i>k */
 401
 402     if (value_notzero_p(pivot)) {
 403       for (i=0; i<rank; i++) {
 404         value_assign(x,H->p[i][k]);
 405         if (value_notzero_p(x)) {
 406           value_modulus(aux,x,pivot);
 407
 408           /* floor[integer division] (corrected for neg x) */
 409           if (value_neg_p(x) && value_notzero_p(aux)) {
 410             value_division(x,x,pivot);
 411             value_decrement(x,x);
 412
 413             /* x=(x/pivot)-1; */
 414           }
 415           else
 416             value_division(x,x,pivot);
 417
 418           /* H->p[i][j] -= x * H->p[rank][j]; */
 419           for (j=0; j<nc; j++) {
 420             value_multiply(aux,x,H->p[rank][j]);
 421             value_subtract(H->p[i][j],H->p[i][j],aux);
 422           }
 423
 424           /* U->p[i][j] -= x * U->p[rank][j]; */
 425           if (U)
 426             for (j=0; j<nr; j++) {
 427               value_multiply(aux,x,U->p[rank][j]);
 428               value_subtract(U->p[i][j],U->p[i][j],aux);
 429             }
 430
 431           /* Q->p[rank][j] += x * Q->p[i][j]; */
 432           if (Q)
 433             for (j=0; j<nr; j++) {
 434               value_addmul(Q->p[rank][j], x, Q->p[i][j]);
 435             }
 436 #ifdef DEBUG
 437           fprintf(stderr,
 438                   "row %d = row %d - %d row %d -----------\n", i+1, i+1, x, rank+1);
 439           Matrix_Print(stderr,0,H);
 440 #endif
 441         } /* if (x) */
 442       } /* for (i) */
 443       rank++;
 444     } /* if (pivot!=0) */
 445   } /* for (k) */
 446
 447   /* Clear all the 'Value' variables */
 448   value_clear(pivot); value_clear(x);
 449   value_clear(aux);
 450   for(i=0;i<nc;i++)
 451     value_clear(temp1[i]);
 452   for(i=0;i<nr;i++)
 453     value_clear(temp2[i]);
 454   free(temp2);
 455   free(temp1);
 456   return rank;
 457 } /* Hermite */
 458
 459 void right_hermite(Matrix *A,Matrix **Hp,Matrix **Up,Matrix **Qp) {
 460
 461   Matrix *H, *Q, *U;
 462   int i, j, nr, nc;
 463   Value tmp;
 464
 465   /* Computes form: A = QH , UA = H */
 466   nc = A->NbColumns;
 467   nr = A->NbRows;
 468
 469   /* H = A */
 470   *Hp = H = Matrix_Alloc(nr,nc);
 471   if (!H) {
 472     errormsg1("DomRightHermite", "outofmem", "out of memory space");
 473     return;
 474   }
 475
 476   /* Initialize all the 'Value' variables */
 477   value_init(tmp);
 478
 479   Vector_Copy(A->p_Init,H->p_Init,nr*nc);
 480
 481   /* U = I */
 482   if (Up) {
 483     *Up = U = Matrix_Alloc(nr, nr);
 484     if (!U) {
 485       errormsg1("DomRightHermite", "outofmem", "out of memory space");
 486       value_clear(tmp);
 487       return;
 488     }
 489     Vector_Set(U->p_Init,0,nr*nr);             /* zero's */
 490     for(i=0;i<nr;i++)                          /* with diagonal of 1's */
 491       value_set_si(U->p[i][i],1);
 492   }
 493   else
 494     U = (Matrix *)0;
 495
 496   /* Q = I */
 497   /* Actually I compute Q transpose... its easier */
 498   if (Qp) {
 499     *Qp = Q = Matrix_Alloc(nr,nr);
 500     if (!Q) {
 501       errormsg1("DomRightHermite", "outofmem", "out of memory space");
 502       value_clear(tmp);
 503       return;
 504     }
 505     Vector_Set(Q->p_Init,0,nr*nr);            /* zero's */
 506     for (i=0;i<nr;i++)                      /* with diagonal of 1's */
 507       value_set_si(Q->p[i][i],1);
 508   }
 509   else
 510     Q = (Matrix *)0;
 511
 512   hermite(H,U,Q);
 513
 514   /* Q is returned transposed */
 515   /* Transpose Q */
 516   if (Q) {
 517     for (i=0; i<nr; i++) {
 518       for (j=i+1; j<nr; j++) {
 519         value_assign(tmp,Q->p[i][j]);
 520         value_assign(Q->p[i][j],Q->p[j][i] );
 521         value_assign(Q->p[j][i],tmp);
 522       }
 523     }
 524   }
 525   value_clear(tmp);
 526   return;
 527 } /* right_hermite */
 528
 529 void left_hermite(Matrix *A,Matrix **Hp,Matrix **Qp,Matrix **Up) {
 530
 531   Matrix *H, *HT, *Q, *U;
 532   int i, j, nc, nr;
 533   Value tmp;
 534
 535   /* Computes left form: A = HQ , AU = H ,
 536                         T    T T    T T   T
 537      using right form  A  = Q H  , U A = H */
 538
 539   nr = A->NbRows;
 540   nc = A->NbColumns;
 541
 542   /* HT = A transpose */
 543   HT = Matrix_Alloc(nc, nr);
 544   if (!HT) {
 545     errormsg1("DomLeftHermite", "outofmem", "out of memory space");
 546     return;
 547   }
 548   value_init(tmp);
 549   for (i=0; i<nr; i++)
 550     for (j=0; j<nc; j++)
 551       value_assign(HT->p[j][i],A->p[i][j]);
 552
 553   /* U = I */
 554   if (Up) {
 555     *Up = U = Matrix_Alloc(nc,nc);
 556     if (!U) {
 557       errormsg1("DomLeftHermite", "outofmem", "out of memory space");
 558       value_clear(tmp);
 559       return;
 560     }
 561     Vector_Set(U->p_Init,0,nc*nc);            /* zero's */
 562     for (i=0;i<nc;i++)                        /* with diagonal of 1's */
 563       value_set_si(U->p[i][i],1);
 564   }
 565   else U=(Matrix *)0;
 566
 567   /* Q = I */
 568   if (Qp) {
 569     *Qp = Q = Matrix_Alloc(nc, nc);
 570     if (!Q) {
 571       errormsg1("DomLeftHermite", "outofmem", "out of memory space");
 572       value_clear(tmp);
 573       return;
 574     }
 575     Vector_Set(Q->p_Init,0,nc*nc);            /* zero's */
 576     for (i=0;i<nc;i++)                        /* with diagonal of 1's */
 577       value_set_si(Q->p[i][i],1);
 578   }
 579   else Q=(Matrix *)0;
 580   hermite(HT,U,Q);
 581
 582   /* H = HT transpose */
 583   *Hp = H = Matrix_Alloc(nr,nc);
 584   if (!H) {
 585     errormsg1("DomLeftHermite", "outofmem", "out of memory space");
 586     value_clear(tmp);
 587     return;
 588   }
 589   for (i=0; i<nr; i++)
 590     for (j=0;j<nc;j++)
 591       value_assign(H->p[i][j],HT->p[j][i]);
 592   Matrix_Free(HT);
 593
 594   /* Transpose U */
 595   if (U) {
 596     for (i=0; i<nc; i++) {
 597       for (j=i+1; j<nc; j++) {
 598         value_assign(tmp,U->p[i][j]);
 599         value_assign(U->p[i][j],U->p[j][i] );
 600         value_assign(U->p[j][i],tmp);
 601       }
 602     }
 603   }
 604   value_clear(tmp);
 605 } /* left_hermite */
 606
 607 /*
 608  * Given a integer matrix 'Mat'(k x k), compute its inverse rational matrix
 609  * 'MatInv' k x (k+1). The last column of each row in matrix MatInv is used
 610  * to store the common denominator of the entries in a row. The output is 1,
 611  * if 'Mat' is non-singular (invertible), otherwise the output is 0. Note::
 612  * (1) Matrix 'Mat' is modified during the inverse operation.
 613  * (2) Matrix 'MatInv' must be preallocated before passing into this function.
 614  */
 615 int MatInverse(Matrix *Mat,Matrix *MatInv ) {
 616
 617   int i, k, j, c;
 618   Value x, gcd, piv;
 619   Value m1,m2;
 620
 621   if(Mat->NbRows != Mat->NbColumns) {
 622    fprintf(stderr,"Trying to invert a non-square matrix !\n");
 623     return 0;
 624   }
 625
 626   /* Initialize all the 'Value' variables */
 627   value_init(x);  value_init(gcd); value_init(piv);
 628   value_init(m1); value_init(m2);
 629
 630   k = Mat->NbRows;
 631
 632   /* Initialise MatInv */
 633   Vector_Set(MatInv->p[0],0,k*(k+1));
 634
 635   /* Initialize 'MatInv' to Identity matrix form. Each diagonal entry is set*/
 636   /* to 1. Last column of each row (denominator of each entry in a row) is  */
 637   /* also set to 1.                                                         */
 638   for(i=0;i<k;++i) {
 639     value_set_si(MatInv->p[i][i],1);
 640     value_set_si(MatInv->p[i][k],1);    /* denum */
 641   }
 642   /* Apply Gauss-Jordan elimination method on the two matrices 'Mat' and  */
 643   /* 'MatInv' in parallel.                                                */
 644   for(i=0;i<k;++i) {
 645
 646     /* Check if the diagonal entry (new pivot) is non-zero or not */
 647     if(value_zero_p(Mat->p[i][i])) {
 648
 649       /* Search for a non-zero pivot down the column(i) */
 650       for(j=i;j<k;++j)
 651         if(value_notzero_p(Mat->p[j][i]))
 652           break;
 653
 654       /* If no non-zero pivot is found, the matrix 'Mat' is non-invertible */
 655       /* Return 0.                                                         */
 656       if(j==k) {
 657
 658         /* Clear all the 'Value' variables */
 659         value_clear(x);  value_clear(gcd); value_clear(piv);
 660         value_clear(m1); value_clear(m2);
 661         return 0;
 662       }
 663
 664       /* Exchange the rows, row(i) and row(j) so that the diagonal element */
 665       /* Mat->p[i][i] (pivot) is non-zero. Repeat the same operations on    */
 666       /* matrix 'MatInv'.                                                   */
 667       for(c=0;c<k;++c) {
 668
 669         /* Interchange rows, row(i) and row(j) of matrix 'Mat'    */
 670         value_assign(x,Mat->p[j][c]);
 671         value_assign(Mat->p[j][c],Mat->p[i][c]);
 672         value_assign(Mat->p[i][c],x);
 673
 674         /* Interchange rows, row(i) and row(j) of matrix 'MatInv' */
 675         value_assign(x,MatInv->p[j][c]);
 676         value_assign(MatInv->p[j][c],MatInv->p[i][c]);
 677         value_assign(MatInv->p[i][c],x);
 678       }
 679     }
 680
 681     /* Make all the entries in column(i) of matrix 'Mat' zero except the */
 682     /* diagonal entry. Repeat the same sequence of operations on matrix  */
 683     /* 'MatInv'.                                                         */
 684     for(j=0;j<k;++j) {
 685       if (j==i) continue;                /* Skip the pivot */
 686       value_assign(x,Mat->p[j][i]);
 687       if(value_notzero_p(x)) {
 688         value_assign(piv,Mat->p[i][i]);
 689         value_gcd(gcd, x, piv);
 690         if (value_notone_p(gcd) ) {
 691           value_divexact(x, x, gcd);
 692           value_divexact(piv, piv, gcd);
 693         }
 694         for(c=((j>i)?i:0);c<k;++c) {
 695           value_multiply(m1,piv,Mat->p[j][c]);
 696           value_multiply(m2,x,Mat->p[i][c]);
 697           value_subtract(Mat->p[j][c],m1,m2);
 698         }
 699         for(c=0;c<k;++c) {
 700           value_multiply(m1,piv,MatInv->p[j][c]);
 701           value_multiply(m2,x,MatInv->p[i][c]);
 702           value_subtract(MatInv->p[j][c],m1,m2);
 703         }
 704
 705         /* Simplify row(j) of the two matrices 'Mat' and 'MatInv' by */
 706         /* dividing the rows with the common GCD.                     */
 707         Vector_Gcd(&MatInv->p[j][0],k,&m1);
 708         Vector_Gcd(&Mat->p[j][0],k,&m2);
 709         value_gcd(gcd, m1, m2);
 710         if(value_notone_p(gcd)) {
 711           for(c=0;c<k;++c) {
 712             value_divexact(Mat->p[j][c], Mat->p[j][c], gcd);
 713             value_divexact(MatInv->p[j][c], MatInv->p[j][c], gcd);
 714           }
 715         }
 716       }
 717     }
 718   }
 719
 720   /* Simplify every row so that 'Mat' reduces to Identity matrix. Perform  */
 721   /* the same sequence of operations on the matrix 'MatInv'.               */
 722   for(j=0;j<k;++j) {
 723     value_assign(MatInv->p[j][k],Mat->p[j][j]);
 724
 725     /* Make the last column (denominator of each entry) of every row greater */
 726     /* than zero.                                                            */
 727     Vector_Normalize_Positive(&MatInv->p[j][0],k+1,k);
 728   }
 729
 730   /* Clear all the 'Value' variables */
 731   value_clear(x);  value_clear(gcd); value_clear(piv);
 732   value_clear(m1); value_clear(m2);
 733
 734   return 1;
 735 } /* Mat_Inverse */
 736
 737 /*
 738  * Given (m x n) integer matrix 'X' and n x (k+1) rational matrix 'P', compute
 739  * the rational m x (k+1) rational matrix  'S'. The last column in each row of
 740  * the rational matrices is used to store the common denominator of elements
 741  * in a row.
 742  */
 743 void rat_prodmat(Matrix *S,Matrix *X,Matrix *P) {
 744
 745   int i,j,k;
 746   int last_column_index = P->NbColumns - 1;
 747   Value lcm, gcd,last_column_entry,s1;
 748   Value m1,m2;
 749
 750   /* Initialize all the 'Value' variables */
 751   value_init(lcm); value_init(gcd);
 752   value_init(last_column_entry); value_init(s1);
 753   value_init(m1); value_init(m2);
 754
 755   /* Compute the LCM of last column entries (denominators) of rows */
 756   value_assign(lcm,P->p[0][last_column_index]);
 757   for(k=1;k<P->NbRows;++k) {
 758     value_assign(last_column_entry,P->p[k][last_column_index]);
 759     value_gcd(gcd, lcm, last_column_entry);
 760     value_divexact(m1, last_column_entry, gcd);
 761     value_multiply(lcm,lcm,m1);
 762   }
 763
 764   /* S[i][j] = Sum(X[i][k] * P[k][j] where Sum is extended over k = 1..nbrows*/
 765   for(i=0;i<X->NbRows;++i)
 766     for(j=0;j<P->NbColumns-1;++j) {
 767
 768       /* Initialize s1 to zero. */
 769       value_set_si(s1,0);
 770       for(k=0;k<P->NbRows;++k) {
 771
 772         /* If the LCM of last column entries is one, simply add the products */
 773         if(value_one_p(lcm)) {
 774           value_addmul(s1, X->p[i][k], P->p[k][j]);
 775         }
 776
 777         /* Numerator (num) and denominator (denom) of S[i][j] is given by :- */
 778         /* numerator  = Sum(X[i][k]*P[k][j]*lcm/P[k][last_column_index]) and */
 779         /* denominator= lcm where Sum is extended over k = 1..nbrows.        */
 780         else {
 781           value_multiply(m1,X->p[i][k],P->p[k][j]);
 782           value_division(m2,lcm,P->p[k][last_column_index]);
 783           value_addmul(s1, m1, m2);
 784         }
 785       }
 786       value_assign(S->p[i][j],s1);
 787     }
 788
 789   for(i=0;i<S->NbRows;++i) {
 790     value_assign(S->p[i][last_column_index],lcm);
 791
 792     /* Normalize the rows so that last element >=0 */
 793     Vector_Normalize_Positive(&S->p[i][0],S->NbColumns,S->NbColumns-1);
 794   }
 795
 796   /* Clear all the 'Value' variables */
 797   value_clear(lcm); value_clear(gcd);
 798   value_clear(last_column_entry); value_clear(s1);
 799   value_clear(m1); value_clear(m2);
 800
 801   return;
 802 } /* rat_prodmat */
 803
 804 /*
 805  * Given a matrix 'Mat' and vector 'p1', compute the matrix-vector product
 806  * and store the result in vector 'p2'.
 807  */
 808 void Matrix_Vector_Product(Matrix *Mat,Value *p1,Value *p2) {
 809
 810   int NbRows, NbColumns, i, j;
 811   Value **cm, *q, *cp1, *cp2;
 812
 813   NbRows=Mat->NbRows;
 814   NbColumns=Mat->NbColumns;
 815
 816   cm = Mat->p;
 817   cp2 = p2;
 818   for(i=0;i<NbRows;i++) {
 819     q = *cm++;
 820     cp1 = p1;
 821     value_multiply(*cp2,*q,*cp1);
 822     q++;
 823     cp1++;
 824
 825     /* *cp2 = *q++ * *cp1++ */
 826     for(j=1;j<NbColumns;j++) {
 827       value_addmul(*cp2, *q, *cp1);
 828       q++;
 829       cp1++;
 830     }
 831     cp2++;
 832   }
 833   return;
 834 } /* Matrix_Vector_Product */
 835
 836 /*
 837  * Given a vector 'p1' and a matrix 'Mat', compute the vector-matrix product
 838  * and store the result in vector 'p2'
 839  */
 840 void Vector_Matrix_Product(Value *p1,Matrix *Mat,Value *p2) {
 841
 842   int NbRows, NbColumns, i, j;
 843   Value **cm, *cp1, *cp2;
 844
 845   NbRows=Mat->NbRows;
 846   NbColumns=Mat->NbColumns;
 847   cp2 = p2;
 848   cm  = Mat->p;
 849   for (j=0;j<NbColumns;j++) {
 850     cp1 = p1;
 851     value_multiply(*cp2,*(*cm+j),*cp1);
 852     cp1++;
 853
 854     /* *cp2= *(*cm+j) * *cp1++; */
 855     for (i=1;i<NbRows;i++) {
 856       value_addmul(*cp2, *(*(cm+i)+j), *cp1);
 857       cp1++;
 858     }
 859     cp2++;
 860   }
 861   return;
 862 } /* Vector_Matrix_Product */
 863
 864 /*
 865  * Given matrices 'Mat1' and 'Mat2', compute the matrix product and store in
 866  * matrix 'Mat3'
 867  */
 868 void Matrix_Product(Matrix *Mat1,Matrix *Mat2,Matrix *Mat3) {
 869
 870   int Size, i, j, k;
 871   unsigned NbRows, NbColumns;
 872   Value **q1, **q2, *p1, *p3,sum;
 873
 874   NbRows    = Mat1->NbRows;
 875   NbColumns = Mat2->NbColumns;
 876
 877   Size      = Mat1->NbColumns;
 878   if(Mat2->NbRows!=Size||Mat3->NbRows!=NbRows||Mat3->NbColumns!=NbColumns) {
 879     fprintf(stderr, "? Matrix_Product : incompatible matrix dimension\n");
 880     return;
 881   }
 882   value_init(sum);
 883   p3 = Mat3->p_Init;
 884   q1 = Mat1->p;
 885   q2 = Mat2->p;
 886
 887   /* Mat3[i][j] = Sum(Mat1[i][k]*Mat2[k][j] where sum is over k = 1..nbrows */
 888   for (i=0;i<NbRows;i++) {
 889     for (j=0;j<NbColumns;j++) {
 890       p1 = *(q1+i);
 891       value_set_si(sum,0);
 892       for (k=0;k<Size;k++) {
 893         value_addmul(sum, *p1, *(*(q2+k)+j));
 894         p1++;
 895       }
 896       value_assign(*p3,sum);
 897       p3++;
 898     }
 899   }
 900   value_clear(sum);
 901   return;
 902 } /* Matrix_Product */
 903
 904 /*
 905  * Given a rational matrix 'Mat'(k x k), compute its inverse rational matrix
 906  * 'MatInv' k x k.
 907  * The output is 1,
 908  * if 'Mat' is non-singular (invertible), otherwise the output is 0. Note::
 909  * (1) Matrix 'Mat' is modified during the inverse operation.
 910  * (2) Matrix 'MatInv' must be preallocated before passing into this function.
 911  */
 912 int Matrix_Inverse(Matrix *Mat,Matrix *MatInv ) {
 913
 914   int i, k, j, c;
 915   Value x, gcd, piv;
 916   Value m1,m2;
 917   Value *den;
 918
 919   if(Mat->NbRows != Mat->NbColumns) {
 920    fprintf(stderr,"Trying to invert a non-square matrix !\n");
 921     return 0;
 922   }
 923
 924   /* Initialize all the 'Value' variables */
 925   value_init(x);  value_init(gcd); value_init(piv);
 926   value_init(m1); value_init(m2);
 927
 928   k = Mat->NbRows;
 929
 930   /* Initialise MatInv */
 931   Vector_Set(MatInv->p[0],0,k*k);
 932
 933   /* Initialize 'MatInv' to Identity matrix form. Each diagonal entry is set*/
 934   /* to 1. Last column of each row (denominator of each entry in a row) is  */
 935   /* also set to 1.                                                         */
 936   for(i=0;i<k;++i) {
 937     value_set_si(MatInv->p[i][i],1);
 938     /* value_set_si(MatInv->p[i][k],1); // denum */
 939   }
 940   /* Apply Gauss-Jordan elimination method on the two matrices 'Mat' and  */
 941   /* 'MatInv' in parallel.                                                */
 942   for(i=0;i<k;++i) {
 943
 944     /* Check if the diagonal entry (new pivot) is non-zero or not */
 945     if(value_zero_p(Mat->p[i][i])) {
 946
 947       /* Search for a non-zero pivot down the column(i) */
 948       for(j=i;j<k;++j)
 949         if(value_notzero_p(Mat->p[j][i]))
 950           break;
 951
 952       /* If no non-zero pivot is found, the matrix 'Mat' is non-invertible */
 953       /* Return 0.                                                         */
 954       if(j==k) {
 955
 956         /* Clear all the 'Value' variables */
 957         value_clear(x);  value_clear(gcd); value_clear(piv);
 958         value_clear(m1); value_clear(m2);
 959         return 0;
 960       }
 961
 962       /* Exchange the rows, row(i) and row(j) so that the diagonal element */
 963       /* Mat->p[i][i] (pivot) is non-zero. Repeat the same operations on    */
 964       /* matrix 'MatInv'.                                                   */
 965       for(c=0;c<k;++c) {
 966
 967         /* Interchange rows, row(i) and row(j) of matrix 'Mat'    */
 968         value_assign(x,Mat->p[j][c]);
 969         value_assign(Mat->p[j][c],Mat->p[i][c]);
 970         value_assign(Mat->p[i][c],x);
 971
 972         /* Interchange rows, row(i) and row(j) of matrix 'MatInv' */
 973         value_assign(x,MatInv->p[j][c]);
 974         value_assign(MatInv->p[j][c],MatInv->p[i][c]);
 975         value_assign(MatInv->p[i][c],x);
 976       }
 977     }
 978
 979     /* Make all the entries in column(i) of matrix 'Mat' zero except the */
 980     /* diagonal entry. Repeat the same sequence of operations on matrix  */
 981     /* 'MatInv'.                                                         */
 982     for(j=0;j<k;++j) {
 983       if (j==i) continue;                /* Skip the pivot */
 984       value_assign(x,Mat->p[j][i]);
 985       if(value_notzero_p(x)) {
 986         value_assign(piv,Mat->p[i][i]);
 987         value_gcd(gcd, x, piv);
 988         if (value_notone_p(gcd) ) {
 989           value_divexact(x, x, gcd);
 990           value_divexact(piv, piv, gcd);
 991         }
 992         for(c=((j>i)?i:0);c<k;++c) {
 993           value_multiply(m1,piv,Mat->p[j][c]);
 994           value_multiply(m2,x,Mat->p[i][c]);
 995           value_subtract(Mat->p[j][c],m1,m2);
 996         }
 997         for(c=0;c<k;++c) {
 998           value_multiply(m1,piv,MatInv->p[j][c]);
 999           value_multiply(m2,x,MatInv->p[i][c]);
1000           value_subtract(MatInv->p[j][c],m1,m2);
1001         }
1002
1003         /* Simplify row(j) of the two matrices 'Mat' and 'MatInv' by */
1004         /* dividing the rows with the common GCD.                     */
1005         Vector_Gcd(&MatInv->p[j][0],k,&m1);
1006         Vector_Gcd(&Mat->p[j][0],k,&m2);
1007         value_gcd(gcd, m1, m2);
1008         if(value_notone_p(gcd)) {
1009           for(c=0;c<k;++c) {
1010             value_divexact(Mat->p[j][c], Mat->p[j][c], gcd);
1011             value_divexact(MatInv->p[j][c], MatInv->p[j][c], gcd);
1012           }
1013         }
1014       }
1015     }
1016   }
1017
1018   /* Find common denom for each row */
1019    den = (Value *)malloc(k*sizeof(Value));
1020    value_set_si(x,1);
1021    for(j=0 ; j<k ; ++j) {
1022      value_init(den[j]);
1023      value_assign(den[j],Mat->p[j][j]);
1024
1025      /* gcd is always positive */
1026      Vector_Gcd(&MatInv->p[j][0],k,&gcd);
1027      value_gcd(gcd, gcd, den[j]);
1028      if (value_neg_p(den[j]))
1029        value_oppose(gcd,gcd); /* make denominator positive */
1030      if (value_notone_p(gcd)) {
1031        for (c=0; c<k; c++)
1032          value_divexact(MatInv->p[j][c], MatInv->p[j][c], gcd); /* normalize */
1033        value_divexact(den[j], den[j], gcd);
1034      }
1035      value_gcd(gcd, x, den[j]);
1036      value_divexact(m1, den[j], gcd);
1037      value_multiply(x,x,m1);
1038    }
1039    if (value_notone_p(x))
1040      for(j=0 ; j<k ; ++j) {
1041        for (c=0; c<k; c++) {
1042          value_division(m1,x,den[j]);
1043          value_multiply(MatInv->p[j][c],MatInv->p[j][c],m1);  /* normalize */
1044        }
1045      }
1046
1047    /* Clear all the 'Value' variables */
1048    for(j=0 ; j<k ; ++j) {
1049      value_clear(den[j]);
1050    }
1051    value_clear(x);  value_clear(gcd); value_clear(piv);
1052    value_clear(m1); value_clear(m2);
1053    free(den);
1054
1055    return 1;
1056 } /* Matrix_Inverse */
1057
1058
1059
1060
1061
1062
1063
1064