Return to Ihermite.c CVS log

Up to [local] / OpenXM / src / Ti

Annotation of OpenXM/src/Ti/Ihermite.c, Revision 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: