]>
Commit | Line | Data |
---|---|---|
98975653 DB |
1 | /* Integer matrix math routines |
2 | Copyright (C) 2003, 2004 Free Software Foundation, Inc. | |
3 | Contributed by Daniel Berlin <dberlin@dberlin.org>. | |
4 | ||
5 | This file is part of GCC. | |
6 | ||
7 | GCC is free software; you can redistribute it and/or modify it under | |
8 | the terms of the GNU General Public License as published by the Free | |
9 | Software Foundation; either version 2, or (at your option) any later | |
10 | version. | |
11 | ||
12 | GCC is distributed in the hope that it will be useful, but WITHOUT ANY | |
13 | WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 | for more details. | |
16 | ||
17 | You should have received a copy of the GNU General Public License | |
18 | along with GCC; see the file COPYING. If not, write to the Free | |
19 | Software Foundation, 59 Temple Place - Suite 330, Boston, MA | |
20 | 02111-1307, USA. */ | |
21 | #include "config.h" | |
22 | #include "system.h" | |
23 | #include "coretypes.h" | |
24 | #include "tm.h" | |
25 | #include "ggc.h" | |
26 | #include "varray.h" | |
36d59cf7 | 27 | #include "tree.h" |
98975653 DB |
28 | #include "lambda.h" |
29 | ||
30 | static void lambda_matrix_get_column (lambda_matrix, int, int, | |
31 | lambda_vector); | |
32 | ||
33 | /* Allocate a matrix of M rows x N cols. */ | |
34 | ||
35 | lambda_matrix | |
36 | lambda_matrix_new (int m, int n) | |
37 | { | |
38 | lambda_matrix mat; | |
39 | int i; | |
40 | ||
41 | mat = ggc_alloc (m * sizeof (lambda_vector)); | |
42 | ||
43 | for (i = 0; i < m; i++) | |
44 | mat[i] = lambda_vector_new (n); | |
45 | ||
46 | return mat; | |
47 | } | |
48 | ||
49 | /* Copy the elements of M x N matrix MAT1 to MAT2. */ | |
50 | ||
51 | void | |
52 | lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2, | |
53 | int m, int n) | |
54 | { | |
55 | int i; | |
56 | ||
57 | for (i = 0; i < m; i++) | |
58 | lambda_vector_copy (mat1[i], mat2[i], n); | |
59 | } | |
60 | ||
61 | /* Store the N x N identity matrix in MAT. */ | |
62 | ||
63 | void | |
64 | lambda_matrix_id (lambda_matrix mat, int size) | |
65 | { | |
66 | int i, j; | |
67 | ||
68 | for (i = 0; i < size; i++) | |
69 | for (j = 0; j < size; j++) | |
70 | mat[i][j] = (i == j) ? 1 : 0; | |
71 | } | |
72 | ||
f67d92e9 DB |
73 | /* Return true if MAT is the identity matrix of SIZE */ |
74 | ||
75 | bool | |
76 | lambda_matrix_id_p (lambda_matrix mat, int size) | |
77 | { | |
78 | int i, j; | |
79 | for (i = 0; i < size; i++) | |
80 | for (j = 0; j < size; j++) | |
81 | { | |
82 | if (i == j) | |
83 | { | |
84 | if (mat[i][j] != 1) | |
85 | return false; | |
86 | } | |
87 | else | |
88 | { | |
89 | if (mat[i][j] != 0) | |
90 | return false; | |
91 | } | |
92 | } | |
93 | return true; | |
94 | } | |
95 | ||
98975653 DB |
96 | /* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */ |
97 | ||
98 | void | |
99 | lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n) | |
100 | { | |
101 | int i; | |
102 | ||
103 | for (i = 0; i < m; i++) | |
104 | lambda_vector_negate (mat1[i], mat2[i], n); | |
105 | } | |
106 | ||
107 | /* Take the transpose of matrix MAT1 and store it in MAT2. | |
108 | MAT1 is an M x N matrix, so MAT2 must be N x M. */ | |
109 | ||
110 | void | |
111 | lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n) | |
112 | { | |
113 | int i, j; | |
114 | ||
115 | for (i = 0; i < n; i++) | |
116 | for (j = 0; j < m; j++) | |
117 | mat2[i][j] = mat1[j][i]; | |
118 | } | |
119 | ||
120 | ||
121 | /* Add two M x N matrices together: MAT3 = MAT1+MAT2. */ | |
122 | ||
123 | void | |
124 | lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2, | |
125 | lambda_matrix mat3, int m, int n) | |
126 | { | |
127 | int i; | |
128 | ||
129 | for (i = 0; i < m; i++) | |
130 | lambda_vector_add (mat1[i], mat2[i], mat3[i], n); | |
131 | } | |
132 | ||
133 | /* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */ | |
134 | ||
135 | void | |
136 | lambda_matrix_add_mc (lambda_matrix mat1, int const1, | |
137 | lambda_matrix mat2, int const2, | |
138 | lambda_matrix mat3, int m, int n) | |
139 | { | |
140 | int i; | |
141 | ||
142 | for (i = 0; i < m; i++) | |
143 | lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n); | |
144 | } | |
145 | ||
146 | /* Multiply two matrices: MAT3 = MAT1 * MAT2. | |
147 | MAT1 is an M x R matrix, and MAT2 is R x N. The resulting MAT2 | |
148 | must therefore be M x N. */ | |
149 | ||
150 | void | |
151 | lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2, | |
152 | lambda_matrix mat3, int m, int r, int n) | |
153 | { | |
154 | ||
155 | int i, j, k; | |
156 | ||
157 | for (i = 0; i < m; i++) | |
158 | { | |
159 | for (j = 0; j < n; j++) | |
160 | { | |
161 | mat3[i][j] = 0; | |
162 | for (k = 0; k < r; k++) | |
163 | mat3[i][j] += mat1[i][k] * mat2[k][j]; | |
164 | } | |
165 | } | |
166 | } | |
167 | ||
168 | /* Get column COL from the matrix MAT and store it in VEC. MAT has | |
169 | N rows, so the length of VEC must be N. */ | |
170 | ||
171 | static void | |
172 | lambda_matrix_get_column (lambda_matrix mat, int n, int col, | |
173 | lambda_vector vec) | |
174 | { | |
175 | int i; | |
176 | ||
177 | for (i = 0; i < n; i++) | |
178 | vec[i] = mat[i][col]; | |
179 | } | |
180 | ||
8c27b7d4 | 181 | /* Delete rows r1 to r2 (not including r2). */ |
98975653 DB |
182 | |
183 | void | |
184 | lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to) | |
185 | { | |
186 | int i; | |
187 | int dist; | |
188 | dist = to - from; | |
189 | ||
190 | for (i = to; i < rows; i++) | |
191 | mat[i - dist] = mat[i]; | |
192 | ||
193 | for (i = rows - dist; i < rows; i++) | |
194 | mat[i] = NULL; | |
195 | } | |
196 | ||
197 | /* Swap rows R1 and R2 in matrix MAT. */ | |
198 | ||
199 | void | |
200 | lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2) | |
201 | { | |
202 | lambda_vector row; | |
203 | ||
204 | row = mat[r1]; | |
205 | mat[r1] = mat[r2]; | |
206 | mat[r2] = row; | |
207 | } | |
208 | ||
209 | /* Add a multiple of row R1 of matrix MAT with N columns to row R2: | |
210 | R2 = R2 + CONST1 * R1. */ | |
211 | ||
212 | void | |
213 | lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1) | |
214 | { | |
215 | int i; | |
216 | ||
217 | if (const1 == 0) | |
218 | return; | |
219 | ||
220 | for (i = 0; i < n; i++) | |
221 | mat[r2][i] += const1 * mat[r1][i]; | |
222 | } | |
223 | ||
224 | /* Negate row R1 of matrix MAT which has N columns. */ | |
225 | ||
226 | void | |
227 | lambda_matrix_row_negate (lambda_matrix mat, int n, int r1) | |
228 | { | |
229 | lambda_vector_negate (mat[r1], mat[r1], n); | |
230 | } | |
231 | ||
232 | /* Multiply row R1 of matrix MAT with N columns by CONST1. */ | |
233 | ||
234 | void | |
235 | lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1) | |
236 | { | |
237 | int i; | |
238 | ||
239 | for (i = 0; i < n; i++) | |
240 | mat[r1][i] *= const1; | |
241 | } | |
242 | ||
243 | /* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */ | |
244 | ||
245 | void | |
246 | lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2) | |
247 | { | |
248 | int i; | |
249 | int tmp; | |
250 | for (i = 0; i < m; i++) | |
251 | { | |
252 | tmp = mat[i][col1]; | |
253 | mat[i][col1] = mat[i][col2]; | |
254 | mat[i][col2] = tmp; | |
255 | } | |
256 | } | |
257 | ||
258 | /* Add a multiple of column C1 of matrix MAT with M rows to column C2: | |
259 | C2 = C2 + CONST1 * C1. */ | |
260 | ||
261 | void | |
262 | lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1) | |
263 | { | |
264 | int i; | |
265 | ||
266 | if (const1 == 0) | |
267 | return; | |
268 | ||
269 | for (i = 0; i < m; i++) | |
270 | mat[i][c2] += const1 * mat[i][c1]; | |
271 | } | |
272 | ||
273 | /* Negate column C1 of matrix MAT which has M rows. */ | |
274 | ||
275 | void | |
276 | lambda_matrix_col_negate (lambda_matrix mat, int m, int c1) | |
277 | { | |
278 | int i; | |
279 | ||
280 | for (i = 0; i < m; i++) | |
281 | mat[i][c1] *= -1; | |
282 | } | |
283 | ||
284 | /* Multiply column C1 of matrix MAT with M rows by CONST1. */ | |
285 | ||
286 | void | |
287 | lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1) | |
288 | { | |
289 | int i; | |
290 | ||
291 | for (i = 0; i < m; i++) | |
292 | mat[i][c1] *= const1; | |
293 | } | |
294 | ||
295 | /* Compute the inverse of the N x N matrix MAT and store it in INV. | |
296 | ||
297 | We don't _really_ compute the inverse of MAT. Instead we compute | |
298 | det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function | |
299 | result. This is necessary to preserve accuracy, because we are dealing | |
300 | with integer matrices here. | |
301 | ||
302 | The algorithm used here is a column based Gauss-Jordan elimination on MAT | |
303 | and the identity matrix in parallel. The inverse is the result of applying | |
304 | the same operations on the identity matrix that reduce MAT to the identity | |
305 | matrix. | |
306 | ||
307 | When MAT is a 2 x 2 matrix, we don't go through the whole process, because | |
308 | it is easily inverted by inspection and it is a very common case. */ | |
309 | ||
310 | static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int); | |
311 | ||
312 | int | |
313 | lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n) | |
314 | { | |
315 | if (n == 2) | |
316 | { | |
317 | int a, b, c, d, det; | |
318 | a = mat[0][0]; | |
319 | b = mat[1][0]; | |
320 | c = mat[0][1]; | |
321 | d = mat[1][1]; | |
322 | inv[0][0] = d; | |
323 | inv[0][1] = -c; | |
324 | inv[1][0] = -b; | |
325 | inv[1][1] = a; | |
326 | det = (a * d - b * c); | |
327 | if (det < 0) | |
328 | { | |
329 | det *= -1; | |
330 | inv[0][0] *= -1; | |
331 | inv[1][0] *= -1; | |
332 | inv[0][1] *= -1; | |
333 | inv[1][1] *= -1; | |
334 | } | |
335 | return det; | |
336 | } | |
337 | else | |
338 | return lambda_matrix_inverse_hard (mat, inv, n); | |
339 | } | |
340 | ||
341 | /* If MAT is not a special case, invert it the hard way. */ | |
342 | ||
343 | static int | |
344 | lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n) | |
345 | { | |
346 | lambda_vector row; | |
347 | lambda_matrix temp; | |
348 | int i, j; | |
349 | int determinant; | |
350 | ||
351 | temp = lambda_matrix_new (n, n); | |
352 | lambda_matrix_copy (mat, temp, n, n); | |
353 | lambda_matrix_id (inv, n); | |
354 | ||
355 | /* Reduce TEMP to a lower triangular form, applying the same operations on | |
356 | INV which starts as the identity matrix. N is the number of rows and | |
357 | columns. */ | |
358 | for (j = 0; j < n; j++) | |
359 | { | |
360 | row = temp[j]; | |
361 | ||
362 | /* Make every element in the current row positive. */ | |
363 | for (i = j; i < n; i++) | |
364 | if (row[i] < 0) | |
365 | { | |
366 | lambda_matrix_col_negate (temp, n, i); | |
367 | lambda_matrix_col_negate (inv, n, i); | |
368 | } | |
369 | ||
370 | /* Sweep the upper triangle. Stop when only the diagonal element in the | |
371 | current row is nonzero. */ | |
372 | while (lambda_vector_first_nz (row, n, j + 1) < n) | |
373 | { | |
374 | int min_col = lambda_vector_min_nz (row, n, j); | |
375 | lambda_matrix_col_exchange (temp, n, j, min_col); | |
376 | lambda_matrix_col_exchange (inv, n, j, min_col); | |
377 | ||
378 | for (i = j + 1; i < n; i++) | |
379 | { | |
380 | int factor; | |
381 | ||
382 | factor = -1 * row[i]; | |
383 | if (row[j] != 1) | |
384 | factor /= row[j]; | |
385 | ||
386 | lambda_matrix_col_add (temp, n, j, i, factor); | |
387 | lambda_matrix_col_add (inv, n, j, i, factor); | |
388 | } | |
389 | } | |
390 | } | |
391 | ||
392 | /* Reduce TEMP from a lower triangular to the identity matrix. Also compute | |
393 | the determinant, which now is simply the product of the elements on the | |
394 | diagonal of TEMP. If one of these elements is 0, the matrix has 0 as an | |
395 | eigenvalue so it is singular and hence not invertible. */ | |
396 | determinant = 1; | |
397 | for (j = n - 1; j >= 0; j--) | |
398 | { | |
399 | int diagonal; | |
400 | ||
401 | row = temp[j]; | |
402 | diagonal = row[j]; | |
403 | ||
41806d92 NS |
404 | /* The matrix must not be singular. */ |
405 | gcc_assert (diagonal); | |
98975653 DB |
406 | |
407 | determinant = determinant * diagonal; | |
408 | ||
409 | /* If the diagonal is not 1, then multiply the each row by the | |
410 | diagonal so that the middle number is now 1, rather than a | |
411 | rational. */ | |
412 | if (diagonal != 1) | |
413 | { | |
414 | for (i = 0; i < j; i++) | |
415 | lambda_matrix_col_mc (inv, n, i, diagonal); | |
416 | for (i = j + 1; i < n; i++) | |
417 | lambda_matrix_col_mc (inv, n, i, diagonal); | |
418 | ||
419 | row[j] = diagonal = 1; | |
420 | } | |
421 | ||
422 | /* Sweep the lower triangle column wise. */ | |
423 | for (i = j - 1; i >= 0; i--) | |
424 | { | |
425 | if (row[i]) | |
426 | { | |
427 | int factor = -row[i]; | |
428 | lambda_matrix_col_add (temp, n, j, i, factor); | |
429 | lambda_matrix_col_add (inv, n, j, i, factor); | |
430 | } | |
431 | ||
432 | } | |
433 | } | |
434 | ||
435 | return determinant; | |
436 | } | |
437 | ||
438 | /* Decompose a N x N matrix MAT to a product of a lower triangular H | |
439 | and a unimodular U matrix such that MAT = H.U. N is the size of | |
440 | the rows of MAT. */ | |
441 | ||
442 | void | |
443 | lambda_matrix_hermite (lambda_matrix mat, int n, | |
444 | lambda_matrix H, lambda_matrix U) | |
445 | { | |
446 | lambda_vector row; | |
447 | int i, j, factor, minimum_col; | |
448 | ||
449 | lambda_matrix_copy (mat, H, n, n); | |
450 | lambda_matrix_id (U, n); | |
451 | ||
452 | for (j = 0; j < n; j++) | |
453 | { | |
454 | row = H[j]; | |
455 | ||
456 | /* Make every element of H[j][j..n] positive. */ | |
457 | for (i = j; i < n; i++) | |
458 | { | |
459 | if (row[i] < 0) | |
460 | { | |
461 | lambda_matrix_col_negate (H, n, i); | |
462 | lambda_vector_negate (U[i], U[i], n); | |
463 | } | |
464 | } | |
465 | ||
8e3c61c5 | 466 | /* Stop when only the diagonal element is nonzero. */ |
98975653 DB |
467 | while (lambda_vector_first_nz (row, n, j + 1) < n) |
468 | { | |
469 | minimum_col = lambda_vector_min_nz (row, n, j); | |
470 | lambda_matrix_col_exchange (H, n, j, minimum_col); | |
471 | lambda_matrix_row_exchange (U, j, minimum_col); | |
472 | ||
473 | for (i = j + 1; i < n; i++) | |
474 | { | |
475 | factor = row[i] / row[j]; | |
476 | lambda_matrix_col_add (H, n, j, i, -1 * factor); | |
477 | lambda_matrix_row_add (U, n, i, j, factor); | |
478 | } | |
479 | } | |
480 | } | |
481 | } | |
482 | ||
483 | /* Given an M x N integer matrix A, this function determines an M x | |
484 | M unimodular matrix U, and an M x N echelon matrix S such that | |
485 | "U.A = S". This decomposition is also known as "right Hermite". | |
486 | ||
487 | Ref: Algorithm 2.1 page 33 in "Loop Transformations for | |
8c27b7d4 | 488 | Restructuring Compilers" Utpal Banerjee. */ |
98975653 DB |
489 | |
490 | void | |
491 | lambda_matrix_right_hermite (lambda_matrix A, int m, int n, | |
492 | lambda_matrix S, lambda_matrix U) | |
493 | { | |
494 | int i, j, i0 = 0; | |
495 | ||
496 | lambda_matrix_copy (A, S, m, n); | |
497 | lambda_matrix_id (U, m); | |
498 | ||
499 | for (j = 0; j < n; j++) | |
500 | { | |
501 | if (lambda_vector_first_nz (S[j], m, i0) < m) | |
502 | { | |
503 | ++i0; | |
504 | for (i = m - 1; i >= i0; i--) | |
505 | { | |
506 | while (S[i][j] != 0) | |
507 | { | |
508 | int sigma, factor, a, b; | |
509 | ||
510 | a = S[i-1][j]; | |
511 | b = S[i][j]; | |
512 | sigma = (a * b < 0) ? -1: 1; | |
513 | a = abs (a); | |
514 | b = abs (b); | |
515 | factor = sigma * (a / b); | |
516 | ||
517 | lambda_matrix_row_add (S, n, i, i-1, -factor); | |
518 | lambda_matrix_row_exchange (S, i, i-1); | |
519 | ||
520 | lambda_matrix_row_add (U, m, i, i-1, -factor); | |
521 | lambda_matrix_row_exchange (U, i, i-1); | |
522 | } | |
523 | } | |
524 | } | |
525 | } | |
526 | } | |
527 | ||
528 | /* Given an M x N integer matrix A, this function determines an M x M | |
529 | unimodular matrix V, and an M x N echelon matrix S such that "A = | |
530 | V.S". This decomposition is also known as "left Hermite". | |
531 | ||
532 | Ref: Algorithm 2.2 page 36 in "Loop Transformations for | |
8c27b7d4 | 533 | Restructuring Compilers" Utpal Banerjee. */ |
98975653 DB |
534 | |
535 | void | |
536 | lambda_matrix_left_hermite (lambda_matrix A, int m, int n, | |
537 | lambda_matrix S, lambda_matrix V) | |
538 | { | |
539 | int i, j, i0 = 0; | |
540 | ||
541 | lambda_matrix_copy (A, S, m, n); | |
542 | lambda_matrix_id (V, m); | |
543 | ||
544 | for (j = 0; j < n; j++) | |
545 | { | |
546 | if (lambda_vector_first_nz (S[j], m, i0) < m) | |
547 | { | |
548 | ++i0; | |
549 | for (i = m - 1; i >= i0; i--) | |
550 | { | |
551 | while (S[i][j] != 0) | |
552 | { | |
553 | int sigma, factor, a, b; | |
554 | ||
555 | a = S[i-1][j]; | |
556 | b = S[i][j]; | |
557 | sigma = (a * b < 0) ? -1: 1; | |
558 | a = abs (a); | |
559 | b = abs (b); | |
560 | factor = sigma * (a / b); | |
561 | ||
562 | lambda_matrix_row_add (S, n, i, i-1, -factor); | |
563 | lambda_matrix_row_exchange (S, i, i-1); | |
564 | ||
565 | lambda_matrix_col_add (V, m, i-1, i, factor); | |
566 | lambda_matrix_col_exchange (V, m, i, i-1); | |
567 | } | |
568 | } | |
569 | } | |
570 | } | |
571 | } | |
572 | ||
8e3c61c5 | 573 | /* When it exists, return the first nonzero row in MAT after row |
98975653 DB |
574 | STARTROW. Otherwise return rowsize. */ |
575 | ||
576 | int | |
577 | lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize, | |
578 | int startrow) | |
579 | { | |
580 | int j; | |
581 | bool found = false; | |
582 | ||
583 | for (j = startrow; (j < rowsize) && !found; j++) | |
584 | { | |
585 | if ((mat[j] != NULL) | |
586 | && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize)) | |
587 | found = true; | |
588 | } | |
589 | ||
590 | if (found) | |
591 | return j - 1; | |
592 | return rowsize; | |
593 | } | |
594 | ||
595 | /* Calculate the projection of E sub k to the null space of B. */ | |
596 | ||
597 | void | |
598 | lambda_matrix_project_to_null (lambda_matrix B, int rowsize, | |
599 | int colsize, int k, lambda_vector x) | |
600 | { | |
601 | lambda_matrix M1, M2, M3, I; | |
602 | int determinant; | |
603 | ||
ea4b7848 | 604 | /* Compute c(I-B^T inv(B B^T) B) e sub k. */ |
98975653 | 605 | |
ea4b7848 | 606 | /* M1 is the transpose of B. */ |
98975653 DB |
607 | M1 = lambda_matrix_new (colsize, colsize); |
608 | lambda_matrix_transpose (B, M1, rowsize, colsize); | |
609 | ||
610 | /* M2 = B * B^T */ | |
611 | M2 = lambda_matrix_new (colsize, colsize); | |
612 | lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize); | |
613 | ||
614 | /* M3 = inv(M2) */ | |
615 | M3 = lambda_matrix_new (colsize, colsize); | |
616 | determinant = lambda_matrix_inverse (M2, M3, rowsize); | |
617 | ||
618 | /* M2 = B^T (inv(B B^T)) */ | |
619 | lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize); | |
620 | ||
621 | /* M1 = B^T (inv(B B^T)) B */ | |
622 | lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize); | |
623 | lambda_matrix_negate (M1, M1, colsize, colsize); | |
624 | ||
625 | I = lambda_matrix_new (colsize, colsize); | |
626 | lambda_matrix_id (I, colsize); | |
627 | ||
628 | lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize); | |
629 | ||
630 | lambda_matrix_get_column (M2, colsize, k - 1, x); | |
631 | ||
632 | } | |
633 | ||
634 | /* Multiply a vector VEC by a matrix MAT. | |
635 | MAT is an M*N matrix, and VEC is a vector with length N. The result | |
636 | is stored in DEST which must be a vector of length M. */ | |
637 | ||
638 | void | |
639 | lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n, | |
640 | lambda_vector vec, lambda_vector dest) | |
641 | { | |
642 | int i, j; | |
643 | ||
644 | lambda_vector_clear (dest, m); | |
645 | for (i = 0; i < m; i++) | |
646 | for (j = 0; j < n; j++) | |
647 | dest[i] += matrix[i][j] * vec[j]; | |
648 | } | |
649 | ||
650 | /* Print out an M x N matrix MAT to OUTFILE. */ | |
651 | ||
652 | void | |
653 | print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n) | |
654 | { | |
655 | int i; | |
656 | ||
657 | for (i = 0; i < m; i++) | |
658 | print_lambda_vector (outfile, matrix[i], n); | |
659 | fprintf (outfile, "\n"); | |
660 | } | |
661 |