1169689Skan/* Integer matrix math routines 2169689Skan Copyright (C) 2003, 2004, 2005 Free Software Foundation, Inc. 3169689Skan Contributed by Daniel Berlin <dberlin@dberlin.org>. 4169689Skan 5169689SkanThis file is part of GCC. 6169689Skan 7169689SkanGCC is free software; you can redistribute it and/or modify it under 8169689Skanthe terms of the GNU General Public License as published by the Free 9169689SkanSoftware Foundation; either version 2, or (at your option) any later 10169689Skanversion. 11169689Skan 12169689SkanGCC is distributed in the hope that it will be useful, but WITHOUT ANY 13169689SkanWARRANTY; without even the implied warranty of MERCHANTABILITY or 14169689SkanFITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 15169689Skanfor more details. 16169689Skan 17169689SkanYou should have received a copy of the GNU General Public License 18169689Skanalong with GCC; see the file COPYING. If not, write to the Free 19169689SkanSoftware Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 20169689Skan02110-1301, USA. */ 21169689Skan#include "config.h" 22169689Skan#include "system.h" 23169689Skan#include "coretypes.h" 24169689Skan#include "tm.h" 25169689Skan#include "ggc.h" 26169689Skan#include "tree.h" 27169689Skan#include "lambda.h" 28169689Skan 29169689Skanstatic void lambda_matrix_get_column (lambda_matrix, int, int, 30169689Skan lambda_vector); 31169689Skan 32169689Skan/* Allocate a matrix of M rows x N cols. */ 33169689Skan 34169689Skanlambda_matrix 35169689Skanlambda_matrix_new (int m, int n) 36169689Skan{ 37169689Skan lambda_matrix mat; 38169689Skan int i; 39169689Skan 40169689Skan mat = ggc_alloc (m * sizeof (lambda_vector)); 41169689Skan 42169689Skan for (i = 0; i < m; i++) 43169689Skan mat[i] = lambda_vector_new (n); 44169689Skan 45169689Skan return mat; 46169689Skan} 47169689Skan 48169689Skan/* Copy the elements of M x N matrix MAT1 to MAT2. */ 49169689Skan 50169689Skanvoid 51169689Skanlambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2, 52169689Skan int m, int n) 53169689Skan{ 54169689Skan int i; 55169689Skan 56169689Skan for (i = 0; i < m; i++) 57169689Skan lambda_vector_copy (mat1[i], mat2[i], n); 58169689Skan} 59169689Skan 60169689Skan/* Store the N x N identity matrix in MAT. */ 61169689Skan 62169689Skanvoid 63169689Skanlambda_matrix_id (lambda_matrix mat, int size) 64169689Skan{ 65169689Skan int i, j; 66169689Skan 67169689Skan for (i = 0; i < size; i++) 68169689Skan for (j = 0; j < size; j++) 69169689Skan mat[i][j] = (i == j) ? 1 : 0; 70169689Skan} 71169689Skan 72169689Skan/* Return true if MAT is the identity matrix of SIZE */ 73169689Skan 74169689Skanbool 75169689Skanlambda_matrix_id_p (lambda_matrix mat, int size) 76169689Skan{ 77169689Skan int i, j; 78169689Skan for (i = 0; i < size; i++) 79169689Skan for (j = 0; j < size; j++) 80169689Skan { 81169689Skan if (i == j) 82169689Skan { 83169689Skan if (mat[i][j] != 1) 84169689Skan return false; 85169689Skan } 86169689Skan else 87169689Skan { 88169689Skan if (mat[i][j] != 0) 89169689Skan return false; 90169689Skan } 91169689Skan } 92169689Skan return true; 93169689Skan} 94169689Skan 95169689Skan/* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */ 96169689Skan 97169689Skanvoid 98169689Skanlambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n) 99169689Skan{ 100169689Skan int i; 101169689Skan 102169689Skan for (i = 0; i < m; i++) 103169689Skan lambda_vector_negate (mat1[i], mat2[i], n); 104169689Skan} 105169689Skan 106169689Skan/* Take the transpose of matrix MAT1 and store it in MAT2. 107169689Skan MAT1 is an M x N matrix, so MAT2 must be N x M. */ 108169689Skan 109169689Skanvoid 110169689Skanlambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n) 111169689Skan{ 112169689Skan int i, j; 113169689Skan 114169689Skan for (i = 0; i < n; i++) 115169689Skan for (j = 0; j < m; j++) 116169689Skan mat2[i][j] = mat1[j][i]; 117169689Skan} 118169689Skan 119169689Skan 120169689Skan/* Add two M x N matrices together: MAT3 = MAT1+MAT2. */ 121169689Skan 122169689Skanvoid 123169689Skanlambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2, 124169689Skan lambda_matrix mat3, int m, int n) 125169689Skan{ 126169689Skan int i; 127169689Skan 128169689Skan for (i = 0; i < m; i++) 129169689Skan lambda_vector_add (mat1[i], mat2[i], mat3[i], n); 130169689Skan} 131169689Skan 132169689Skan/* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */ 133169689Skan 134169689Skanvoid 135169689Skanlambda_matrix_add_mc (lambda_matrix mat1, int const1, 136169689Skan lambda_matrix mat2, int const2, 137169689Skan lambda_matrix mat3, int m, int n) 138169689Skan{ 139169689Skan int i; 140169689Skan 141169689Skan for (i = 0; i < m; i++) 142169689Skan lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n); 143169689Skan} 144169689Skan 145169689Skan/* Multiply two matrices: MAT3 = MAT1 * MAT2. 146169689Skan MAT1 is an M x R matrix, and MAT2 is R x N. The resulting MAT2 147169689Skan must therefore be M x N. */ 148169689Skan 149169689Skanvoid 150169689Skanlambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2, 151169689Skan lambda_matrix mat3, int m, int r, int n) 152169689Skan{ 153169689Skan 154169689Skan int i, j, k; 155169689Skan 156169689Skan for (i = 0; i < m; i++) 157169689Skan { 158169689Skan for (j = 0; j < n; j++) 159169689Skan { 160169689Skan mat3[i][j] = 0; 161169689Skan for (k = 0; k < r; k++) 162169689Skan mat3[i][j] += mat1[i][k] * mat2[k][j]; 163169689Skan } 164169689Skan } 165169689Skan} 166169689Skan 167169689Skan/* Get column COL from the matrix MAT and store it in VEC. MAT has 168169689Skan N rows, so the length of VEC must be N. */ 169169689Skan 170169689Skanstatic void 171169689Skanlambda_matrix_get_column (lambda_matrix mat, int n, int col, 172169689Skan lambda_vector vec) 173169689Skan{ 174169689Skan int i; 175169689Skan 176169689Skan for (i = 0; i < n; i++) 177169689Skan vec[i] = mat[i][col]; 178169689Skan} 179169689Skan 180169689Skan/* Delete rows r1 to r2 (not including r2). */ 181169689Skan 182169689Skanvoid 183169689Skanlambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to) 184169689Skan{ 185169689Skan int i; 186169689Skan int dist; 187169689Skan dist = to - from; 188169689Skan 189169689Skan for (i = to; i < rows; i++) 190169689Skan mat[i - dist] = mat[i]; 191169689Skan 192169689Skan for (i = rows - dist; i < rows; i++) 193169689Skan mat[i] = NULL; 194169689Skan} 195169689Skan 196169689Skan/* Swap rows R1 and R2 in matrix MAT. */ 197169689Skan 198169689Skanvoid 199169689Skanlambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2) 200169689Skan{ 201169689Skan lambda_vector row; 202169689Skan 203169689Skan row = mat[r1]; 204169689Skan mat[r1] = mat[r2]; 205169689Skan mat[r2] = row; 206169689Skan} 207169689Skan 208169689Skan/* Add a multiple of row R1 of matrix MAT with N columns to row R2: 209169689Skan R2 = R2 + CONST1 * R1. */ 210169689Skan 211169689Skanvoid 212169689Skanlambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1) 213169689Skan{ 214169689Skan int i; 215169689Skan 216169689Skan if (const1 == 0) 217169689Skan return; 218169689Skan 219169689Skan for (i = 0; i < n; i++) 220169689Skan mat[r2][i] += const1 * mat[r1][i]; 221169689Skan} 222169689Skan 223169689Skan/* Negate row R1 of matrix MAT which has N columns. */ 224169689Skan 225169689Skanvoid 226169689Skanlambda_matrix_row_negate (lambda_matrix mat, int n, int r1) 227169689Skan{ 228169689Skan lambda_vector_negate (mat[r1], mat[r1], n); 229169689Skan} 230169689Skan 231169689Skan/* Multiply row R1 of matrix MAT with N columns by CONST1. */ 232169689Skan 233169689Skanvoid 234169689Skanlambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1) 235169689Skan{ 236169689Skan int i; 237169689Skan 238169689Skan for (i = 0; i < n; i++) 239169689Skan mat[r1][i] *= const1; 240169689Skan} 241169689Skan 242169689Skan/* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */ 243169689Skan 244169689Skanvoid 245169689Skanlambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2) 246169689Skan{ 247169689Skan int i; 248169689Skan int tmp; 249169689Skan for (i = 0; i < m; i++) 250169689Skan { 251169689Skan tmp = mat[i][col1]; 252169689Skan mat[i][col1] = mat[i][col2]; 253169689Skan mat[i][col2] = tmp; 254169689Skan } 255169689Skan} 256169689Skan 257169689Skan/* Add a multiple of column C1 of matrix MAT with M rows to column C2: 258169689Skan C2 = C2 + CONST1 * C1. */ 259169689Skan 260169689Skanvoid 261169689Skanlambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1) 262169689Skan{ 263169689Skan int i; 264169689Skan 265169689Skan if (const1 == 0) 266169689Skan return; 267169689Skan 268169689Skan for (i = 0; i < m; i++) 269169689Skan mat[i][c2] += const1 * mat[i][c1]; 270169689Skan} 271169689Skan 272169689Skan/* Negate column C1 of matrix MAT which has M rows. */ 273169689Skan 274169689Skanvoid 275169689Skanlambda_matrix_col_negate (lambda_matrix mat, int m, int c1) 276169689Skan{ 277169689Skan int i; 278169689Skan 279169689Skan for (i = 0; i < m; i++) 280169689Skan mat[i][c1] *= -1; 281169689Skan} 282169689Skan 283169689Skan/* Multiply column C1 of matrix MAT with M rows by CONST1. */ 284169689Skan 285169689Skanvoid 286169689Skanlambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1) 287169689Skan{ 288169689Skan int i; 289169689Skan 290169689Skan for (i = 0; i < m; i++) 291169689Skan mat[i][c1] *= const1; 292169689Skan} 293169689Skan 294169689Skan/* Compute the inverse of the N x N matrix MAT and store it in INV. 295169689Skan 296169689Skan We don't _really_ compute the inverse of MAT. Instead we compute 297169689Skan det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function 298169689Skan result. This is necessary to preserve accuracy, because we are dealing 299169689Skan with integer matrices here. 300169689Skan 301169689Skan The algorithm used here is a column based Gauss-Jordan elimination on MAT 302169689Skan and the identity matrix in parallel. The inverse is the result of applying 303169689Skan the same operations on the identity matrix that reduce MAT to the identity 304169689Skan matrix. 305169689Skan 306169689Skan When MAT is a 2 x 2 matrix, we don't go through the whole process, because 307169689Skan it is easily inverted by inspection and it is a very common case. */ 308169689Skan 309169689Skanstatic int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int); 310169689Skan 311169689Skanint 312169689Skanlambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n) 313169689Skan{ 314169689Skan if (n == 2) 315169689Skan { 316169689Skan int a, b, c, d, det; 317169689Skan a = mat[0][0]; 318169689Skan b = mat[1][0]; 319169689Skan c = mat[0][1]; 320169689Skan d = mat[1][1]; 321169689Skan inv[0][0] = d; 322169689Skan inv[0][1] = -c; 323169689Skan inv[1][0] = -b; 324169689Skan inv[1][1] = a; 325169689Skan det = (a * d - b * c); 326169689Skan if (det < 0) 327169689Skan { 328169689Skan det *= -1; 329169689Skan inv[0][0] *= -1; 330169689Skan inv[1][0] *= -1; 331169689Skan inv[0][1] *= -1; 332169689Skan inv[1][1] *= -1; 333169689Skan } 334169689Skan return det; 335169689Skan } 336169689Skan else 337169689Skan return lambda_matrix_inverse_hard (mat, inv, n); 338169689Skan} 339169689Skan 340169689Skan/* If MAT is not a special case, invert it the hard way. */ 341169689Skan 342169689Skanstatic int 343169689Skanlambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n) 344169689Skan{ 345169689Skan lambda_vector row; 346169689Skan lambda_matrix temp; 347169689Skan int i, j; 348169689Skan int determinant; 349169689Skan 350169689Skan temp = lambda_matrix_new (n, n); 351169689Skan lambda_matrix_copy (mat, temp, n, n); 352169689Skan lambda_matrix_id (inv, n); 353169689Skan 354169689Skan /* Reduce TEMP to a lower triangular form, applying the same operations on 355169689Skan INV which starts as the identity matrix. N is the number of rows and 356169689Skan columns. */ 357169689Skan for (j = 0; j < n; j++) 358169689Skan { 359169689Skan row = temp[j]; 360169689Skan 361169689Skan /* Make every element in the current row positive. */ 362169689Skan for (i = j; i < n; i++) 363169689Skan if (row[i] < 0) 364169689Skan { 365169689Skan lambda_matrix_col_negate (temp, n, i); 366169689Skan lambda_matrix_col_negate (inv, n, i); 367169689Skan } 368169689Skan 369169689Skan /* Sweep the upper triangle. Stop when only the diagonal element in the 370169689Skan current row is nonzero. */ 371169689Skan while (lambda_vector_first_nz (row, n, j + 1) < n) 372169689Skan { 373169689Skan int min_col = lambda_vector_min_nz (row, n, j); 374169689Skan lambda_matrix_col_exchange (temp, n, j, min_col); 375169689Skan lambda_matrix_col_exchange (inv, n, j, min_col); 376169689Skan 377169689Skan for (i = j + 1; i < n; i++) 378169689Skan { 379169689Skan int factor; 380169689Skan 381169689Skan factor = -1 * row[i]; 382169689Skan if (row[j] != 1) 383169689Skan factor /= row[j]; 384169689Skan 385169689Skan lambda_matrix_col_add (temp, n, j, i, factor); 386169689Skan lambda_matrix_col_add (inv, n, j, i, factor); 387169689Skan } 388169689Skan } 389169689Skan } 390169689Skan 391169689Skan /* Reduce TEMP from a lower triangular to the identity matrix. Also compute 392169689Skan the determinant, which now is simply the product of the elements on the 393169689Skan diagonal of TEMP. If one of these elements is 0, the matrix has 0 as an 394169689Skan eigenvalue so it is singular and hence not invertible. */ 395169689Skan determinant = 1; 396169689Skan for (j = n - 1; j >= 0; j--) 397169689Skan { 398169689Skan int diagonal; 399169689Skan 400169689Skan row = temp[j]; 401169689Skan diagonal = row[j]; 402169689Skan 403169689Skan /* The matrix must not be singular. */ 404169689Skan gcc_assert (diagonal); 405169689Skan 406169689Skan determinant = determinant * diagonal; 407169689Skan 408169689Skan /* If the diagonal is not 1, then multiply the each row by the 409169689Skan diagonal so that the middle number is now 1, rather than a 410169689Skan rational. */ 411169689Skan if (diagonal != 1) 412169689Skan { 413169689Skan for (i = 0; i < j; i++) 414169689Skan lambda_matrix_col_mc (inv, n, i, diagonal); 415169689Skan for (i = j + 1; i < n; i++) 416169689Skan lambda_matrix_col_mc (inv, n, i, diagonal); 417169689Skan 418169689Skan row[j] = diagonal = 1; 419169689Skan } 420169689Skan 421169689Skan /* Sweep the lower triangle column wise. */ 422169689Skan for (i = j - 1; i >= 0; i--) 423169689Skan { 424169689Skan if (row[i]) 425169689Skan { 426169689Skan int factor = -row[i]; 427169689Skan lambda_matrix_col_add (temp, n, j, i, factor); 428169689Skan lambda_matrix_col_add (inv, n, j, i, factor); 429169689Skan } 430169689Skan 431169689Skan } 432169689Skan } 433169689Skan 434169689Skan return determinant; 435169689Skan} 436169689Skan 437169689Skan/* Decompose a N x N matrix MAT to a product of a lower triangular H 438169689Skan and a unimodular U matrix such that MAT = H.U. N is the size of 439169689Skan the rows of MAT. */ 440169689Skan 441169689Skanvoid 442169689Skanlambda_matrix_hermite (lambda_matrix mat, int n, 443169689Skan lambda_matrix H, lambda_matrix U) 444169689Skan{ 445169689Skan lambda_vector row; 446169689Skan int i, j, factor, minimum_col; 447169689Skan 448169689Skan lambda_matrix_copy (mat, H, n, n); 449169689Skan lambda_matrix_id (U, n); 450169689Skan 451169689Skan for (j = 0; j < n; j++) 452169689Skan { 453169689Skan row = H[j]; 454169689Skan 455169689Skan /* Make every element of H[j][j..n] positive. */ 456169689Skan for (i = j; i < n; i++) 457169689Skan { 458169689Skan if (row[i] < 0) 459169689Skan { 460169689Skan lambda_matrix_col_negate (H, n, i); 461169689Skan lambda_vector_negate (U[i], U[i], n); 462169689Skan } 463169689Skan } 464169689Skan 465169689Skan /* Stop when only the diagonal element is nonzero. */ 466169689Skan while (lambda_vector_first_nz (row, n, j + 1) < n) 467169689Skan { 468169689Skan minimum_col = lambda_vector_min_nz (row, n, j); 469169689Skan lambda_matrix_col_exchange (H, n, j, minimum_col); 470169689Skan lambda_matrix_row_exchange (U, j, minimum_col); 471169689Skan 472169689Skan for (i = j + 1; i < n; i++) 473169689Skan { 474169689Skan factor = row[i] / row[j]; 475169689Skan lambda_matrix_col_add (H, n, j, i, -1 * factor); 476169689Skan lambda_matrix_row_add (U, n, i, j, factor); 477169689Skan } 478169689Skan } 479169689Skan } 480169689Skan} 481169689Skan 482169689Skan/* Given an M x N integer matrix A, this function determines an M x 483169689Skan M unimodular matrix U, and an M x N echelon matrix S such that 484169689Skan "U.A = S". This decomposition is also known as "right Hermite". 485169689Skan 486169689Skan Ref: Algorithm 2.1 page 33 in "Loop Transformations for 487169689Skan Restructuring Compilers" Utpal Banerjee. */ 488169689Skan 489169689Skanvoid 490169689Skanlambda_matrix_right_hermite (lambda_matrix A, int m, int n, 491169689Skan lambda_matrix S, lambda_matrix U) 492169689Skan{ 493169689Skan int i, j, i0 = 0; 494169689Skan 495169689Skan lambda_matrix_copy (A, S, m, n); 496169689Skan lambda_matrix_id (U, m); 497169689Skan 498169689Skan for (j = 0; j < n; j++) 499169689Skan { 500169689Skan if (lambda_vector_first_nz (S[j], m, i0) < m) 501169689Skan { 502169689Skan ++i0; 503169689Skan for (i = m - 1; i >= i0; i--) 504169689Skan { 505169689Skan while (S[i][j] != 0) 506169689Skan { 507169689Skan int sigma, factor, a, b; 508169689Skan 509169689Skan a = S[i-1][j]; 510169689Skan b = S[i][j]; 511169689Skan sigma = (a * b < 0) ? -1: 1; 512169689Skan a = abs (a); 513169689Skan b = abs (b); 514169689Skan factor = sigma * (a / b); 515169689Skan 516169689Skan lambda_matrix_row_add (S, n, i, i-1, -factor); 517169689Skan lambda_matrix_row_exchange (S, i, i-1); 518169689Skan 519169689Skan lambda_matrix_row_add (U, m, i, i-1, -factor); 520169689Skan lambda_matrix_row_exchange (U, i, i-1); 521169689Skan } 522169689Skan } 523169689Skan } 524169689Skan } 525169689Skan} 526169689Skan 527169689Skan/* Given an M x N integer matrix A, this function determines an M x M 528169689Skan unimodular matrix V, and an M x N echelon matrix S such that "A = 529169689Skan V.S". This decomposition is also known as "left Hermite". 530169689Skan 531169689Skan Ref: Algorithm 2.2 page 36 in "Loop Transformations for 532169689Skan Restructuring Compilers" Utpal Banerjee. */ 533169689Skan 534169689Skanvoid 535169689Skanlambda_matrix_left_hermite (lambda_matrix A, int m, int n, 536169689Skan lambda_matrix S, lambda_matrix V) 537169689Skan{ 538169689Skan int i, j, i0 = 0; 539169689Skan 540169689Skan lambda_matrix_copy (A, S, m, n); 541169689Skan lambda_matrix_id (V, m); 542169689Skan 543169689Skan for (j = 0; j < n; j++) 544169689Skan { 545169689Skan if (lambda_vector_first_nz (S[j], m, i0) < m) 546169689Skan { 547169689Skan ++i0; 548169689Skan for (i = m - 1; i >= i0; i--) 549169689Skan { 550169689Skan while (S[i][j] != 0) 551169689Skan { 552169689Skan int sigma, factor, a, b; 553169689Skan 554169689Skan a = S[i-1][j]; 555169689Skan b = S[i][j]; 556169689Skan sigma = (a * b < 0) ? -1: 1; 557169689Skan a = abs (a); 558169689Skan b = abs (b); 559169689Skan factor = sigma * (a / b); 560169689Skan 561169689Skan lambda_matrix_row_add (S, n, i, i-1, -factor); 562169689Skan lambda_matrix_row_exchange (S, i, i-1); 563169689Skan 564169689Skan lambda_matrix_col_add (V, m, i-1, i, factor); 565169689Skan lambda_matrix_col_exchange (V, m, i, i-1); 566169689Skan } 567169689Skan } 568169689Skan } 569169689Skan } 570169689Skan} 571169689Skan 572169689Skan/* When it exists, return the first nonzero row in MAT after row 573169689Skan STARTROW. Otherwise return rowsize. */ 574169689Skan 575169689Skanint 576169689Skanlambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize, 577169689Skan int startrow) 578169689Skan{ 579169689Skan int j; 580169689Skan bool found = false; 581169689Skan 582169689Skan for (j = startrow; (j < rowsize) && !found; j++) 583169689Skan { 584169689Skan if ((mat[j] != NULL) 585169689Skan && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize)) 586169689Skan found = true; 587169689Skan } 588169689Skan 589169689Skan if (found) 590169689Skan return j - 1; 591169689Skan return rowsize; 592169689Skan} 593169689Skan 594169689Skan/* Calculate the projection of E sub k to the null space of B. */ 595169689Skan 596169689Skanvoid 597169689Skanlambda_matrix_project_to_null (lambda_matrix B, int rowsize, 598169689Skan int colsize, int k, lambda_vector x) 599169689Skan{ 600169689Skan lambda_matrix M1, M2, M3, I; 601169689Skan int determinant; 602169689Skan 603169689Skan /* Compute c(I-B^T inv(B B^T) B) e sub k. */ 604169689Skan 605169689Skan /* M1 is the transpose of B. */ 606169689Skan M1 = lambda_matrix_new (colsize, colsize); 607169689Skan lambda_matrix_transpose (B, M1, rowsize, colsize); 608169689Skan 609169689Skan /* M2 = B * B^T */ 610169689Skan M2 = lambda_matrix_new (colsize, colsize); 611169689Skan lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize); 612169689Skan 613169689Skan /* M3 = inv(M2) */ 614169689Skan M3 = lambda_matrix_new (colsize, colsize); 615169689Skan determinant = lambda_matrix_inverse (M2, M3, rowsize); 616169689Skan 617169689Skan /* M2 = B^T (inv(B B^T)) */ 618169689Skan lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize); 619169689Skan 620169689Skan /* M1 = B^T (inv(B B^T)) B */ 621169689Skan lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize); 622169689Skan lambda_matrix_negate (M1, M1, colsize, colsize); 623169689Skan 624169689Skan I = lambda_matrix_new (colsize, colsize); 625169689Skan lambda_matrix_id (I, colsize); 626169689Skan 627169689Skan lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize); 628169689Skan 629169689Skan lambda_matrix_get_column (M2, colsize, k - 1, x); 630169689Skan 631169689Skan} 632169689Skan 633169689Skan/* Multiply a vector VEC by a matrix MAT. 634169689Skan MAT is an M*N matrix, and VEC is a vector with length N. The result 635169689Skan is stored in DEST which must be a vector of length M. */ 636169689Skan 637169689Skanvoid 638169689Skanlambda_matrix_vector_mult (lambda_matrix matrix, int m, int n, 639169689Skan lambda_vector vec, lambda_vector dest) 640169689Skan{ 641169689Skan int i, j; 642169689Skan 643169689Skan lambda_vector_clear (dest, m); 644169689Skan for (i = 0; i < m; i++) 645169689Skan for (j = 0; j < n; j++) 646169689Skan dest[i] += matrix[i][j] * vec[j]; 647169689Skan} 648169689Skan 649169689Skan/* Print out an M x N matrix MAT to OUTFILE. */ 650169689Skan 651169689Skanvoid 652169689Skanprint_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n) 653169689Skan{ 654169689Skan int i; 655169689Skan 656169689Skan for (i = 0; i < m; i++) 657169689Skan print_lambda_vector (outfile, matrix[i], n); 658169689Skan fprintf (outfile, "\n"); 659169689Skan} 660169689Skan 661