1)abbrev package UGAUSS UnitGaussianElimination
2++ References:
3++   J-G. Dumas, F. Heckenbach, D. Sauders, V. Welker:
4++   Computing simplicial homology based on efficient Smith
5++   normal form algorithms.
6++ Description: Gaussian elimination using only 1 or -1 as pivots
7++
8
9UnitGaussianElimination : with
10    pre_gauss : Matrix(Integer) -> Matrix(Integer)
11      ++ pre_gauss(m) performs Gaussian eliminaton on rows of m
12      ++ using only 1 and -1 as pivots.
13      ++ Note:  m is modified in place
14    pre_smith : Matrix(Integer) -> Matrix(Integer)
15      ++ pre_smith(m) first performs pre_gauss(m) and then
16      ++ drops all rows and colums of m contaning pivots.
17    pre_lr : (Matrix(Integer), Vector(Integer), Vector(Integer)
18             ) -> Matrix(Integer)
19      ++ pre_lr(m, pi, pj) is like pre_gauss, but records positions
20      ++ of pivots in pi and pj and returns transformation matrix.
21      ++ m, pi, pj are modified in place
22  == add
23
24    M ==> Matrix(Integer)
25    V ==> Vector(Integer)
26    NNI ==> NonNegativeInteger
27
28    -- Gaussian elimination using only 1 or -1 as pivots
29    -- m is input matrix,
30    -- pivotsj is position of pivot within row, 0 if none
31    -- pivotsk is position of pivot within column, 0 if none
32    -- pivotsj and pivotsk are 0 on entry and filled with
33    -- appropriate values on exit
34    pre_lr0(m, pivotsj, pivotsk, adjust) ==>
35        j := nrows(m)
36        k := ncols(m)
37        for i in 1..j repeat
38            -- choose pivot in row i
39            for l in 1..k repeat
40                pivotsk(l) ~= 0 => "iterate"
41                pv := qelt(m, i, l)
42                if pv = 1 or pv = -1 then
43                     pivotsj(i) := l
44                     pivotsk(l) := i
45                     break
46            pivotsj(i) = 0 => "iterate"
47            l0 := pivotsj(i)
48            ml := qelt(m, i, l0)
49            -- zero out column l0, do not bother with rows containing pivots
50            for n in 1..j repeat
51                if qelt(pivotsj, n) = 0 and qelt(m, n, l0) ~= 0 then
52                    pp := -ml*qelt(m, n, l0)
53                    for l in 1..k repeat
54                        qsetelt!(m, n, l, qelt(m, n, l) + pp*qelt(m, i, l))
55                    adjust
56
57    pre_gauss0(m : M, pivotsj : V, pivotsk : V) : Void ==
58        pre_lr0(m, pivotsj, pivotsk, 0)
59
60    pre_lr(m : M, pivotsj : V, pivotsk : V) : M ==
61        m2 := scalarMatrix(nrows(m), 1)$M
62        adjust ==>
63            for l in 1..j repeat
64                qsetelt!(m2, n, l, qelt(m2, n, l) + pp*qelt(m2, i, l))
65        pre_lr0(m, pivotsj, pivotsk, adjust)
66        m2
67
68    pre_gauss(m : M) : M ==
69        j := nrows(m)
70        k := ncols(m)
71        pivotsj := zero(j)$V
72        pivotsk := zero(k)$V
73        pre_gauss0(m, pivotsj, pivotsk)
74        m
75
76    pre_smith(m : M) : M ==
77        j := nrows(m)
78        k := ncols(m)
79        pivotsj := zero(j)$V
80        pivotsk := zero(k)$V
81        pre_gauss0(m, pivotsj, pivotsk)
82
83        -- count rows containing pivots
84        count := 0
85        for i in 1..j repeat
86            if pivotsj(i) ~= 0 then
87                count := count + 1
88        count = 0 => m
89        -- skip rows and columns containing pivots
90        ress := zero((j - count)::NNI, (k - count)::NNI)$M
91        i0 := 0
92        for i in 1..j repeat
93            if pivotsj(i) = 0 then
94                i0 := i0 + 1
95                l0 := 0
96                for l in 1..k repeat
97                    if qelt(pivotsk, l) = 0 then
98                        l0 := l0 + 1
99                        qsetelt!(ress, i0, l0, qelt(m, i, l))
100        ress
101
102
103)abbrev package ISMITH IntegerSmithNormalForm
104
105M ==> Matrix(Integer)
106SmithForm ==> Record(Smith : M, leftEqMat : M, rightEqMat : M)
107IntegerSmithNormalForm : with
108
109    smith : M -> M
110      ++ smith(m) computes Smith normal form of m
111    smith : (M, M -> M) -> M
112      ++ smith(m, full) computes Smith normal form of m.  It
113      ++ first reduces m to smaller matrix and then uses full
114      ++ to finish.
115    completeSmith : (M, M -> SmithForm) -> SmithForm
116      ++ completeSmith(m, full) computes record containing Smith normal
117      ++ form of m and the left and right equivalence matrices.
118      ++ It first reduces m to smaller matrix and then uses full
119      ++ to finish.
120
121  == add
122
123    PS ==> UnitGaussianElimination
124    V ==> Vector(Integer)
125    SP ==> SmithNormalForm(Integer, V, V, M)
126
127    smith(m : M, full : M -> M) : M ==
128        m1 := copy(m)
129        m2 := pre_smith(m1)$PS
130        m3 := full(m2)
131        j := nrows(m)
132        k := ncols(m)
133        nn := min(j, k)
134        count := j - nrows(m2)
135        res := zero(j, k)$M
136        for i in 1..count repeat
137            qsetelt!(res, i, i, 1)
138        nn2 := nn - count
139        for i in 1..nn2 repeat
140            qsetelt!(res, i + count, i + count, qelt(m3, i, i))
141        res
142
143    smith(m : M) : M == smith(m, smith$SP)
144
145    completeSmith(m : M, full : M -> SmithForm) : SmithForm ==
146        j := nrows(m)
147        k := ncols(m)
148        j > k =>
149            res0 := completeSmith(transpose(m), full)
150            [transpose(res0.Smith), transpose(res0.rightEqMat),
151              transpose(res0.leftEqMat)]
152        pivotsj := zero(j)$V
153        pivotsk := zero(k)$V
154        m1 := copy(m)
155        m2 := pre_lr(m1, pivotsj, pivotsk)$PS
156        lj0 : List(Integer) := [i for i in 1..j | pivotsj(i) = 0]
157        lj1 : List(Integer) := [s for i in 1..j | (s := pivotsj(i)) ~= 0]
158        lj2 : List(Integer) := [i for i in 1..j | pivotsj(i) ~= 0]
159        lk0 : List(Integer) := [i for i in 1..k | pivotsk(i) = 0]
160        lk : List(Integer) := [s for i in 1..k | (s := pivotsk(i)) ~= 0]
161        ljs := concat(lj2, lj0)
162        lks := concat(lj1, lk0)
163        m3 := m1(ljs, lks)
164        m4 := m2(ljs, 1..j)
165        count := #lj2
166        m5 := m3((count + 1)..j, (count + 1)..k)
167        res1 := full(m5)
168        m6 := m4((count + 1)..j, 1..j)
169        m7 := res1.leftEqMat*m6
170        m8 := vertConcat(m4(1..count, 1..j), m7)
171        m9 := m3((count + 1)..j, 1..k)
172        m10 := res1.leftEqMat*m9
173        m11 := vertConcat(m3(1..count, 1..k), m10)
174        m12 := zero(k, k)$M
175        for i in 1..k for l in lks repeat
176            m12(l, i) := 1
177        m13 := m12(1..k, (count + 1)..k)*res1.rightEqMat
178        m14 := horizConcat(m12(1..k, 1..count), m13)
179        m15 := m11(1..j, (count + 1)..k)*res1.rightEqMat
180        m16 := horizConcat(m11(1..j, 1..count), m15)
181        for i in count..1 by -1 repeat
182            if qelt(m16, i, i) = -1 then
183                for l in 1..i repeat
184                    qsetelt!(m16, l, i, -qelt(m16, l, i))
185                for l in 1..k repeat
186                    qsetelt!(m14, l, i, -qelt(m14, l, i))
187            qelt(m16, i, i) ~= 1 =>
188                error "completeSmith: wrong diagonal"
189            for l in 1..(i - 1) repeat
190                pp := qelt(m16, l, i)
191                for n in (i + 1)..k repeat
192                    qsetelt!(m16, l, n, qelt(m16, l, n) - pp*qelt(m16, i, n))
193            for l in 1..k repeat
194                pp := qelt(m14, l, i)
195                for n in (i + 1)..k repeat
196                    qsetelt!(m14, l, n, qelt(m14, l, n) - pp*qelt(m16, i, n))
197            for n in (i + 1)..k repeat qsetelt!(m16, i, n,  0)
198        [m16, m8, m14]
199