Return to Ihermite.c CVS log

Up to [local] / OpenXM / src / Ti

Annotation of OpenXM/src/Ti/Ihermite.c, Revision 1.1.1.1

1.1       maekawa     1: /*
                      2: ** Ihermite.c                                 Birk Huber, 2/99
                      3: **   -- compute hermite normal form of an integer matrix
                      4: **
                      5: **
                      6: ** TiGERS,  Toric Groebner Basis Enumeration by Reverse Search
                      7: ** copyright (c) 1999  Birk Huber
                      8: **
                      9: */
                     10: #include<stdio.h>
                     11: #include<stdlib.h>
                     12: #include "utils.h"
                     13:
                     14:
                     15: /*
                     16: **
                     17: ** Imatrix_hermite:
                     18: **    Input: S an n by m Imatrix.
                     19: **   Output: True
                     20: **
                     21: **   Sidefects: U points to a unitary nxn matrix
                     22: **              S gets replaced by its hermite normal form (U*S)
                     23: **
                     24: **  Method:
                     25: **    U = Inxn
                     26: **    c = 1
                     27: **    while (c<=n) do
                     28: **      find mc>=c such that abs(S(mc,c)) has minimum non-zero value
                     29: **      if (mc!=c) then  switch rows mc and c (in S, and U)
                     30: **      if (S(c,c)<0) then multiply row c by -1
                     31: **      if (S(c,c)!=0) then
                     32: **        for i in c+1:n  add S(i,c)/S(c,c) times row c to row i;
                     33: **      end(if)
                     34: **      if S(i,c)==0 for all i in c+1 ... n set c=c+1
                     35: **    end(while)
                     36: **
                     37: */
                     38: #define Mref(U,i,j) U[(i)-1][(j)-1]
                     39: #define Mrows(U) m
                     40: #define Mcols(U) n
                     41:
                     42:
                     43:
                     44: int ihermite(int **S,int **U,int m, int n){
                     45:   int c=1,mc,i,j,done,sign;
                     46:   int t=0,mv=0, mn,crk=0;
                     47:
                     48:   if (m>n) mn=n; else mn=m;
                     49:
                     50:   /* Initialize U to nxn identity */
                     51:   for(i=1;i<=n;i++){
                     52:     for(j=1;j<=n;j++){
                     53:       if (i==j) Mref(U,i,j)=1;
                     54:       else Mref(U,i,j)=0;
                     55:     }
                     56:   }
                     57:
                     58:   while(c<=mn){
                     59:     /* find minimum entry in col c */
                     60:     mv=Mref(S,c,c);
                     61:     if (mv<0) mv*=-1;
                     62:     mc=c;
                     63:     for(i=c+1;i<=n;i++){
                     64:       t=Mref(S,i,c);
                     65:       if (t<0) t*=-1;
                     66:       if(mv==0 || (mv>t && t!=0)){
                     67:            mv=t;
                     68:            mc=i;
                     69:       }
                     70:     }
                     71:
                     72:     /* if nescesary pivot to put min in row c and multiply by+-1
                     73:        to ensure diagonal entry is positive */
                     74:     if (mc!=c||Mref(S,mc,c)<0){
                     75:       if (Mref(S,mc,c)<0) sign=-1;
                     76:       else sign=+1;
                     77:       for(j=c;j<=m;j++){
                     78:         t=Mref(S,mc,j);
                     79:         Mref(S,mc,j)=Mref(S,c,j);
                     80:         Mref(S,c,j)=sign*t;
                     81:       }
                     82:       for(j=1;j<=n;j++){
                     83:         t=Mref(U,mc,j);
                     84:         Mref(U,mc,j)=Mref(U,c,j);
                     85:         Mref(U,c,j)=sign*t;
                     86:       }
                     87:     }
                     88:
                     89:     /* if collumb is not zero do a reduction step */
                     90:     done=TRUE;
                     91:     if (Mref(S,c,c)!=0){
                     92:       crk=c;
                     93:       for(i=c+1;i<=n;i++){
                     94:         t=Mref(S,i,c)/Mref(S,c,c);
                     95:         for(j=c;j<=m;j++){
                     96:            Mref(S,i,j)-=t*Mref(S,c,j);
                     97:         }
                     98:         for(j=1;j<=n;j++){
                     99:            Mref(U,i,j)-=t*Mref(U,c,j);
                    100:         }
                    101:         if (Mref(S,i,c)!=0) done=FALSE;
                    102:       }
                    103:     }
                    104:     /* if all entrees of col c bellow row c are zero go tonext col */
                    105:     if (done==TRUE) c++;
                    106:   }
                    107:   return crk;
                    108:  }
                    109:
                    110:
                    111:
                    112:
                    113:
                    114:
                    115:
                    116: