source: git/kernel/linear_algebra/linearAlgebra.h @ 23a78e

spielwiese
Last change on this file since 23a78e was 23a78e, checked in by Adi Popescu <adi_popescum@…>, 10 years ago
Separating Headers: kernel/linear_algebra
  • Property mode set to 100644
File size: 25.9 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file lineareAlgebra.h
5 *
6 * This file provides basic linear algebra functionality.
7 *
8 * ABSTRACT: This file provides basic algorithms from linear algebra over
9 * any SINGULAR-supported field.
10 * For the time being, the procedures defined in this file expect matrices
11 * containing objects of the SINGULAR type 'number'. This means that, when
12 * 'currentRing' represents some polynomial ring K[x_1, x_2, ..., x_n], then
13 * the entries of the matrices are 'numbers' representing elements of K (and
14 * NOT 'polys' in K[x_1, x_2, ..., x_n]).
15 * This restriction may become obselete in the future.
16 *
17 * @author Frank Seelisch
18 *
19 *
20 **/
21/*****************************************************************************/
22
23#ifndef LINEAR_ALGEBRA_H
24#define LINEAR_ALGEBRA_H
25
26// include basic SINGULAR structures
27#include <kernel/structs.h>
28
29/**
30 * LU-decomposition of a given (m x n)-matrix.
31 *
32 * Given an (m x n) matrix A, the method computes its LU-decomposition,
33 * that is, it computes matrices P, L, and U such that<br>
34 * - P * A = L * U,<br>
35 * - P is an (m x m) permutation matrix, i.e., its row/columns form the
36 *   standard basis of K^m,<br>
37 * - L is an (m x m) matrix in lower triangular form with all diagonal
38 *   entries equal to 1,<br>
39 * - U is an (m x n) matrix in upper row echelon form.<br>
40 * From these conditions, it follows immediately that also A = P * L * U,
41 * since P is self-inverse.<br>
42 *
43 * The method will modify all argument matrices except aMat, so that
44 * they will finally contain the matrices P, L, and U as specified
45 * above.
46 **/
47void luDecomp(
48       const matrix aMat, /**< [in]  the initial matrix A,          */
49       matrix &pMat,      /**< [out] the row permutation matrix P   */
50       matrix &lMat,      /**< [out] the lower triangular matrix L  */
51       matrix &uMat,      /**< [out] the upper row echelon matrix U */
52       const ring r= currRing       /**< [in] current ring */
53             );
54
55/**
56 * Returns a pivot score for any non-zero matrix entry.
57 *
58 * The smaller the score the better will n serve as a pivot element
59 * in subsequent Gauss elimination steps.
60 *
61 * @return the pivot score of n
62 **/
63int pivotScore(
64               number n,    /**< [in] a non-zero matrix entry */
65               const ring r= currRing /**< [in] current ring */
66              );
67
68/**
69 * Finds the best pivot element in some part of a given matrix.
70 *
71 * Given any matrix A with valid row indices r1..r2 and valid column
72 * indices c1..c2, this method finds the best pivot element for
73 * subsequent Gauss elimination steps in A[r1..r2, c1..c2]. "Best"
74 * here means best with respect to the implementation of the method
75 * 'pivotScore(number n)'.<br>
76 * In the case that all elements in A[r1..r2, c1..c2] are zero, the
77 * method returns false, otherwise true.
78 *
79 * @return false if all relevant matrix entries are zero, true otherwise
80 * @sa pivotScore(number n)
81 **/
82bool pivot(
83           const matrix aMat, /**< [in]  any matrix with number entries */
84           const int r1,      /**< [in]  lower row index                */
85           const int r2,      /**< [in]  upper row index                */
86           const int c1,      /**< [in]  lower column index             */
87           const int c2,      /**< [in]  upper column index             */
88           int* bestR,        /**< [out] address of row index of best
89                                         pivot element                  */
90           int* bestC,        /**< [out] address of column index of
91                                         best pivot element             */
92           const ring r= currRing       /**< [in] current ring */
93          );
94
95/**
96 * Computes the inverse of a given (n x n)-matrix.
97 *
98 * This method expects an (n x n)-matrix, that is, it must have as many
99 * rows as columns. Inversion of the first argument is attempted via the
100 * LU-decomposition. There are two cases:<br>
101 * 1) The matrix is not invertible. Then the method returns false, and
102 *    &iMat remains unchanged.<br>
103 * 2) The matrix is invertible. Then the method returns true, and the
104 *    content of iMat is the inverse of aMat.
105 *
106 * @return true iff aMat is invertible, false otherwise
107 * @sa luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
108 *                           const matrix uMat, matrix &iMat)
109 **/
110bool luInverse(
111               const matrix aMat, /**< [in]  matrix to be inverted */
112               matrix &iMat,      /**< [out] inverse of aMat if
113                                             invertible            */
114               const ring r=currRing /**< [in] current ring */
115              );
116
117/**
118 * Computes the inverse of a given (n x n)-matrix in upper right
119 * triangular form.
120 *
121 * This method expects a quadratic matrix, that is, it must have as
122 * many rows as columns. Moreover, uMat[i, j] = 0, at least for all
123 * i > j, that is, u is in upper right triangular form.<br>
124 * If the argument diagonalIsOne is true, then we know additionally,
125 * that uMat[i, i] = 1, for all i. In this case uMat is invertible.
126 * Contrariwise, if diagonalIsOne is false, we do not know anything
127 * about the diagonal entries. (Note that they may still all be
128 * 1.)<br>
129 * In general, there are two cases:<br>
130 * 1) The matrix is not invertible. Then the method returns false,
131 *    and &iMat remains unchanged.<br>
132 * 2) The matrix is invertible. Then the method returns true, and
133 *    the content of iMat is the inverse of uMat.
134 *
135 * @return true iff uMat is invertible, false otherwise
136 **/
137bool upperRightTriangleInverse(
138       const matrix uMat, /**< [in]  (n x n)-matrix in upper right
139                                     triangular form               */
140       matrix &iMat,      /**< [out] inverse of uMat if invertible */
141       bool diagonalIsOne,/**< [in]  if true, then all diagonal
142                                     entries of uMat are 1         */
143       const ring r= currRing /**< [in] current ring */
144                              );
145
146/**
147 * Computes the inverse of a given (n x n)-matrix in lower left
148 * triangular form.
149 *
150 * This method expects an (n x n)-matrix, that is, it must have as
151 * many rows as columns. Moreover, lMat[i,j] = 0, at least for all
152 * j > i, that ism lMat is in lower left triangular form.<br>
153 * If the argument diagonalIsOne is true, then we know additionally,
154 * that lMat[i, i] = 1, for all i. In this case lMat is invertible.
155 * Contrariwise, if diagonalIsOne is false, we do not know anything
156 * about the diagonal entries. (Note that they may still all be
157 * 1.)<br>
158 * In general, there are two cases:<br>
159 * 1) The matrix is not invertible. Then the method returns false,
160 *    and &iMat remains unchanged.<br>
161 * 2) The matrix is invertible. Then the method returns true, and
162 *    the content of iMat is the inverse of lMat.
163 *
164 * @return true iff lMat is invertible, false otherwise
165 **/
166bool lowerLeftTriangleInverse(
167       const matrix lMat, /**< [in]  (n x n)-matrix in lower left
168                                     triangular form               */
169       matrix &iMat,      /**< [out] inverse of lMat if invertible */
170       bool diagonalIsOne /**< [in]  if true, then all diagonal
171                                     entries of lMat are 1         */
172                              );
173
174/**
175 * Computes the inverse of an (n x n)-matrix which is given by its LU-
176 * decomposition.
177 *
178 * With A denoting the matrix to be inverted, the method expects the
179 * LU-decomposition of A, that is, pMat * A = lMat * uMat, where
180 * the argument matrices have the appropriate proteries; see method
181 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
182 * matrix &uMat)'.<br>
183 * Furthermore, uMat is expected to be an (n x n)-matrix. Then A^(-1)
184 * is computed according to A^(-1) = uMat^(-1) * lMat^(-1) * pMat,
185 * since pMat is self-inverse. This will work if and only if uMat is
186 * invertible, because lMat and pMat are in any case invertible.<br>
187 * In general, there are two cases:<br>
188 * 1) uMat and hence A is not invertible. Then the method returns
189 *    false, and &iMat remains unchanged.<br>
190 * 2) uMat and hence A is invertible. Then the method returns true,
191 *    and the content of iMat is the inverse of A.
192 *
193 * @return true if A is invertible, false otherwise
194 * @sa luInverse(const matrix aMat, matrix &iMat)
195 **/
196bool luInverseFromLUDecomp(
197       const matrix pMat, /**< [in]  permutation matrix of an LU-
198                                     decomposition                */
199       const matrix lMat, /**< [in]  lower left matrix of an LU-
200                                     decomposition                */
201       const matrix uMat, /**< [in]  upper right matrix of an LU-
202                                     decomposition                */
203       matrix &iMat,      /**< [out] inverse of A if invertible   */
204       const ring r= currRing
205                          );
206
207/**
208 * Computes the rank of a given (m x n)-matrix.
209 *
210 * The matrix may already be given in row echelon form, e.g., when
211 * the user has before called luDecomp and passes the upper triangular
212 * matrix U to luRank. In this case, the second argument can be set to
213 * true in order to pass this piece of information to the method.
214 * Otherwise, luDecomp will be called first to compute the matrix U.
215 * The rank is then read off the matrix U.
216 *
217 * @return the rank as an int
218 * @sa luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
219 **/
220int luRank(
221       const matrix aMat,      /**< [in]  input matrix */
222       const bool isRowEchelon,/**< [in]  if true then aMat is assumed to be
223                                          already in row echelon form */
224       const ring r= currRing
225          );
226
227/**
228 * Solves the linear system A * x = b, where A is an (m x n)-matrix
229 * which is given by its LU-decomposition.
230 *
231 * The method expects the LU-decomposition of A, that is,
232 * pMat * A = lMat * uMat, where the argument matrices have the
233 * appropriate proteries; see method
234 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
235 * matrix &uMat)'.<br>
236 * Instead of trying to invert A and return x = A^(-1)*b, this
237 * method
238 * 1) computes b' = pMat * b,
239 * 2) solves the simple system L * y = b', and then
240 * 3) solves the simple system U * x = y.
241 * Note that steps 1) and 2) will always work, as L is in any case
242 * invertible. Moreover, y is uniquely determined. Step 3) will only
243 * work if and only if y is in the column span of U. In that case,
244 * there may be infinitely many solutions.
245 * Thus, there are three cases:<br>
246 * 1) There is no solution. Then the method returns false, and &xVec
247 *    as well as &H remain unchanged.<br>
248 * 2) There is a unique solution. Then the method returns true and
249 *    H is the 1x1 matrix with zero entry.<br>
250 * 3) There are infinitely many solutions. Then the method returns
251 *    true and some solution of the given original linear system.
252 *    Furthermore, the columns of H span the solution space of the
253 *    homogeneous linear system. The dimension of the solution
254 *    space is then the number of columns of H.
255 *
256 * @return true if there is at least one solution, false otherwise
257 **/
258bool luSolveViaLUDecomp(
259       const matrix pMat, /**< [in]  permutation matrix of an LU-
260                                     decomposition                */
261       const matrix lMat, /**< [in]  lower left matrix of an LU-
262                                     decomposition                */
263       const matrix uMat, /**< [in]  upper right matrix of an LU-
264                                     decomposition                */
265       const matrix bVec, /**< [in]  right-hand side of the linear
266                                     system to be solved          */
267       matrix &xVec,      /**< [out] solution of A*x = b          */
268       matrix &H          /**< [out] matrix with columns spanning
269                                     homogeneous solution space   */
270                          );
271
272/**
273 * Solves the linear system A * x = b, where A is an (m x n)-matrix
274 * which is given by its LDU-decomposition.
275 *
276 * The method expects the LDU-decomposition of A, that is,
277 * pMat * A = lMat * dMat^(-1) * uMat, where the argument matrices have
278 * the appropriate proteries; see method
279 * 'lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
280 * matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)'.<br>
281 * Instead of trying to invert A and return x = A^(-1)*b, this
282 * method
283 * 1) computes b' = l * pMat * b,
284 * 2) solves the simple system L * y = b',
285 * 3) computes y' = u * dMat * y,
286 * 4) solves the simple system U * y'' = y',
287 * 5) computes x = 1/(lTimesU) * y''.
288 * Note that steps 1), 2) and 3) will always work, as L is in any case
289 * invertible. Moreover, y and thus y' are uniquely determined.
290 * Step 4) will only work if and only if y' is in the column span of U.
291 * In that case, there may be infinitely many solutions.
292 * In contrast to luSolveViaLUDecomp, this algorithm guarantees that
293 * all divisions which have to be performed in steps 2) and 4) will
294 * work without remainder. This is due to properties of the given LDU-
295 * decomposition. Only the final step 5) performs a division of a vector
296 * by a member of the ground field. (Suppose the ground field is Q or
297 * some rational function field. Then, if A contains only integers or
298 * polynomials, respectively, then steps 1) - 4) will keep all data
299 * integer or polynomial, respectively. This may speed up computations
300 * considerably.)
301 * For the outcome, there are three cases:<br>
302 * 1) There is no solution. Then the method returns false, and &xVec
303 *    as well as &H remain unchanged.<br>
304 * 2) There is a unique solution. Then the method returns true and
305 *    H is the 1x1 matrix with zero entry.<br>
306 * 3) There are infinitely many solutions. Then the method returns
307 *    true and some solution of the given original linear system.
308 *    Furthermore, the columns of H span the solution space of the
309 *    homogeneous linear system. The dimension of the solution
310 *    space is then the number of columns of H.
311 *
312 * @return true if there is at least one solution, false otherwise
313 **/
314bool luSolveViaLDUDecomp(
315       const matrix pMat,  /**< [in]  permutation matrix of an LDU-
316                                      decomposition                     */
317       const matrix lMat,  /**< [in]  lower left matrix of an LDU-
318                                      decomposition                     */
319       const matrix dMat,  /**< [in]  diagonal matrix of an LDU-
320                                      decomposition                     */
321       const matrix uMat,  /**< [in]  upper right matrix of an LDU-
322                                      decomposition                     */
323       const poly l,       /**< [in]  pivot product l of an LDU decomp. */
324       const poly u,       /**< [in]  pivot product u of an LDU decomp. */
325       const poly lTimesU, /**< [in]  product of l and u                */
326       const matrix bVec,  /**< [in]  right-hand side of the linear
327                                      system to be solved               */
328       matrix &xVec,       /**< [out] solution of A*x = b               */
329       matrix &H           /**< [out] matrix with columns spanning
330                                      homogeneous solution space        */
331                          );
332
333/**
334 * Creates a new matrix which is the (nxn) unit matrix, and returns true
335 * in case of success.
336 *
337 * The method will be successful whenever n >= 1. In this case and only then
338 * true will be returned and the new (nxn) unit matrix will reside inside
339 * the second argument.
340 *
341 * @return true iff the (nxn) unit matrix could be build
342 **/
343bool unitMatrix(
344       const int n,     /**< [in]  size of the matrix */
345       matrix &unitMat,  /**< [out] the new (nxn) unit matrix */
346       const ring r= currRing /** [in] current ring */
347               );
348
349/**
350 * Creates a new matrix which is a submatrix of the first argument, and
351 * returns true in case of success.
352 *
353 * The method will be successful whenever rowIndex1 <= rowIndex2 and
354 * colIndex1 <= colIndex2. In this case and only then true will be
355 * returned and the last argument will afterwards contain a copy of the
356 * respective submatrix of the first argument.
357 * Note that this method may also be used to copy an entire matrix.
358 *
359 * @return true iff the submatrix could be build
360 **/
361bool subMatrix(
362       const matrix aMat,    /**< [in]  the input matrix */
363       const int rowIndex1,  /**< [in]  lower row index */
364       const int rowIndex2,  /**< [in]  higher row index */
365       const int colIndex1,  /**< [in]  lower column index */
366       const int colIndex2,  /**< [in]  higher column index */
367       matrix &subMat        /**< [out] the new submatrix */
368              );
369
370/**
371 * Swaps two rows of a given matrix and thereby modifies it.
372 *
373 * The method expects two valid, distinct row indices, i.e. in
374 * [1..n], where n is the number of rows in aMat.
375 **/
376void swapRows(
377       int row1,     /**< [in]     index of first row to swap */
378       int row2,     /**< [in]     index of second row to swap */
379       matrix& aMat  /**< [in/out] matrix subject to swapping */
380             );
381
382/**
383 * Swaps two columns of a given matrix and thereby modifies it.
384 *
385 * The method expects two valid, distinct column indices, i.e. in
386 * [1..n], where n is the number of columns in aMat.
387 **/
388void swapColumns(
389       int column1,  /**< [in]     index of first column to swap */
390       int column2,  /**< [in]     index of second column to swap */
391       matrix& aMat  /**< [in/out] matrix subject to swapping */
392             );
393
394/**
395 * Creates a new matrix which contains the first argument in the top left
396 * corner, and the second in the bottom right.
397 *
398 * All other entries of the resulting matrix which will be created in the
399 * third argument, are zero.
400 **/
401void matrixBlock(
402       const matrix aMat,    /**< [in]  the top left input matrix */
403       const matrix bMat,    /**< [in]  the bottom right input matrix */
404       matrix &block         /**< [out] the new block matrix */
405                );
406
407/**
408 * Computes the characteristic polynomial from a quadratic (2x2) matrix
409 * and returns true in case of success.
410 *
411 * The method will be successful whenever the input matrix is a (2x2) matrix.
412 * In this case, the resulting polynomial will be a univariate polynomial in
413 * the first ring variable of degree 2 and it will reside in the second
414 * argument.
415 * The method assumes that the matrix entries are all numbers, i.e., elements
416 * from the ground field/ring.
417 *
418 * @return true iff the input matrix was (2x2)
419 **/
420bool charPoly(
421       const matrix aMat,    /**< [in]  the input matrix */
422       poly &charPoly        /**< [out] the characteristic polynomial */
423             );
424
425/**
426 * Computes the square root of a non-negative real number and returns
427 * it as a new number.
428 *
429 * The method assumes that the current ground field is the complex
430 * numbers, and that the argument has imaginary part zero.
431 * If the real part is negative, then false is returned, otherwise true.
432 * The square root will be computed in the last argument by Heron's
433 * iteration formula with the first argument as the starting value. The
434 * iteration will stop as soon as two successive values (in the resulting
435 * sequence) differ by no more than the given tolerance, which is assumed
436 * to be a positive real number.
437 *
438 * @return the square root of the non-negative argument number
439 **/
440bool realSqrt(
441       const number n,          /**< [in]  the input number */
442       const number tolerance,  /**< [in]  accuracy of iteration */
443       number &root             /**< [out] the root of n */
444             );
445
446/**
447 * Computes the Hessenberg form of a given square matrix.
448 *
449 * The method assumes that all matrix entries are numbers coming from some
450 * subfield of the reals..
451 * Afterwards, the following conditions will hold:
452 * 1) hessenbergMat = pMat * aMat * pMat^{-1},
453 * 2) hessenbergMat is in Hessenberg form.
454 * During the algorithm, pMat will be constructed as the product of self-
455 * inverse matrices.
456 * The algorithm relies on computing square roots of floating point numbers.
457 * These will be combuted by Heron's iteration formula, with iteration
458 * stopping when two successive approximations of the square root differ by
459 * no more than the given tolerance, which is assumed to be a positve real
460 * number.
461 **/
462void hessenberg(
463       const matrix aMat,      /**< [in]  the square input matrix */
464       matrix &pMat,           /**< [out] the transformation matrix */
465       matrix &hessenbergMat,  /**< [out] the Hessenberg form of aMat */
466       const number tolerance, /**< [in]  accuracy */
467       const ring r
468               );
469
470/**
471 * Returns all solutions of a quadratic polynomial relation with real
472 * coefficients.
473 *
474 * The method assumes that the polynomial is univariate in the first
475 * ring variable, and that the current ground field is the complex numbers.
476 * The polynomial has degree <= 2. Thus, there may be
477 * a) infinitely many zeros, when p == 0; then -1 is returned;
478 * b) no zero, when p == constant <> 0; then 0 is returned;
479 * c) one zero, when p is linear; then 1 is returned;
480 * d) one double zero; then 2 is returned;
481 * e) two distinct zeros; then 3 is returned;
482 * In the case e), s1 and s2 will contain the two distinct solutions.
483 * In cases c) and d) s1 will contain the single/double solution.
484 *
485 * @return the number of distinct zeros
486 **/
487int quadraticSolve(
488       const poly p,           /**< [in]  the polynomial        */
489       number &s1,             /**< [out] first zero, if any    */
490       number &s2,             /**< [out] second zero, if any   */
491       const number tolerance  /**< [in] accuracy               */
492                  );
493
494
495/**
496 * Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a
497 * certain degree in x, whenever a factorization of h(0, y) is given.
498 *
499 * The algorithm is based on Hensel's lemma: Let h(x, y) denote a monic
500 * polynomial in y of degree m + n with coefficients in K[[x]]. Suppose there
501 * are two monic factors f_0(y) (of degree n) and g_0(y) of degree (m) such
502 * that h(0, y) = f_0(y) * g_0(y) and <f_0, g_0> = K[y]. Fix an integer d >= 0.
503 * Then there are monic polynomials in y with coefficients in K[[x]], namely
504 * f(x, y) of degree n and g(x, y) of degree m such that
505 *    h(x, y) = f(x, y) * g(x, y) modulo <x^(d+1)>   (*).
506 *
507 * This implementation expects h, f0, g0, and d as described and computes the
508 * factors f and g. Effectively, h will be given as an element of K[x, y] since
509 * all terms of h with x-degree larger than d can be ignored due to (*).
510 * The method expects the ground ring to contain at least two variables; then
511 * x is the first ring variable, specified by the input index xIndex, and y the
512 * second one, specified by yIndex.
513 *
514 * This code was placed here since the algorithm works by successively solving
515 * d linear equation systems. It is hence an application of other methods
516 * defined in this h-file and its corresponding cc-file.
517 *
518 **/
519void henselFactors(
520       const int xIndex,  /**< [in]  index of x in {1, ..., nvars(basering)} */
521       const int yIndex,  /**< [in]  index of y in {1, ..., nvars(basering)} */
522       const poly h,      /**< [in]  the polynomial h(x, y) as above         */
523       const poly f0,     /**< [in]  the first univariate factor of h(0, y)  */
524       const poly g0,     /**< [in]  the second univariate factor of h(0, y) */
525       const int d,       /**< [in]  the degree bound, d >= 0                */
526       poly &f,           /**< [out] the first factor of h(x, y)             */
527       poly &g            /**< [out] the second factor of h(x, y)            */
528                              );
529
530/**
531 * LU-decomposition of a given (m x n)-matrix with performing only those
532 * divisions that yield zero remainders.
533 *
534 * Given an (m x n) matrix A, the method computes its LDU-decomposition,
535 * that is, it computes matrices P, L, D, and U such that<br>
536 * - P * A = L * D^(-1) * U,<br>
537 * - P is an (m x m) permutation matrix, i.e., its row/columns form the
538 *   standard basis of K^m,<br>
539 * - L is an (m x m) matrix in lower triangular form of full rank,<br>
540 * - D is an (m x m) diagonal matrix with full rank, and<br>
541 * - U is an (m x n) matrix in upper row echelon form.<br>
542 * From these conditions, it follows immediately that also
543 * A = P * L * D^(-1) * U, since P is self-inverse.<br>
544 *
545 * In contrast to luDecomp, this method only performs those divisions which
546 * yield zero remainders. Hence, when e.g. computing over a rational function
547 * field and starting with polynomial entries only (or over Q and starting
548 * with integer entries only), then any newly computed matrix entry will again
549 * be polynomial (or integer).
550 *
551 * The method will modify all argument matrices except aMat, so that
552 * they will finally contain the matrices P, L, D, and U as specified
553 * above. Moreover, in order to further speed up subsequent calls of
554 * luSolveViaLDUDecomp, the two denominators and their product will also be
555 * returned.
556 **/
557void lduDecomp(
558       const matrix aMat, /**< [in]  the initial matrix A,          */
559       matrix &pMat,      /**< [out] the row permutation matrix P   */
560       matrix &lMat,      /**< [out] the lower triangular matrix L  */
561       matrix &dMat,      /**< [out] the diagonal matrix D          */
562       matrix &uMat,      /**< [out] the upper row echelon matrix U */
563       poly &l,           /**< [out] the product of pivots of L     */
564       poly &u,           /**< [out] the product of pivots of U     */
565       poly &lTimesU      /**< [out] the product of l and u         */
566             );
567
568/* helper for qrDoubleShift */
569bool qrDS( const int n, matrix* queue, int& queueL, number* eigenValues,
570       int& eigenValuesL, const number tol1, const number tol2, const ring R);
571
572/**
573 * Tries to find the number n in the array nn[0..nnLength-1].
574 **/
575int similar(
576       const number* nn,       /**< [in] array of numbers to look-up */
577       const int nnLength,     /**< [in] length of nn                */
578       const number n,         /**< [in] number to loop-up in nn     */
579       const number tolerance  /**< [in] tolerance for comparison    */
580           );
581#endif
582/* LINEAR_ALGEBRA_H */
Note: See TracBrowser for help on using the repository browser.