My Project
linearAlgebra.h
Go to the documentation of this file.
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  **/
47 void 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  **/
63 int 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  **/
82 bool 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  **/
110 bool 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  **/
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  **/
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  **/
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  **/
220 int 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  **/
258 bool 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  **/
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  **/
343 bool 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  **/
361 bool 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  **/
376 void 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  **/
388 void 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  **/
401 void 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  **/
420 bool 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  **/
440 bool 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 positive real
460  * number.
461  **/
462 void 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  **/
487 int 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  **/
519 void 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  **/
557 void 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 */
569 bool 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  **/
575 int 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 */
int l
Definition: cfEzgcd.cc:100
int p
Definition: cfModGcd.cc:4080
g
Definition: cfModGcd.cc:4092
FILE * f
Definition: checklibs.c:9
CanonicalForm H
Definition: facAbsFact.cc:60
STATIC_VAR Poly * h
Definition: janet.cc:971
bool luSolveViaLDUDecomp(const matrix pMat, const matrix lMat, const matrix dMat, const matrix uMat, const poly l, const poly u, const poly lTimesU, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LDU-decomposit...
bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat, const ring r=currRing)
Computes the inverse of an (n x n)-matrix which is given by its LU- decomposition.
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring r=currRing)
LU-decomposition of a given (m x n)-matrix.
bool realSqrt(const number n, const number tolerance, number &root)
Computes the square root of a non-negative real number and returns it as a new number.
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
bool luInverse(const matrix aMat, matrix &iMat, const ring r=currRing)
Computes the inverse of a given (n x n)-matrix.
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring r)
Computes the Hessenberg form of a given square matrix.
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
bool unitMatrix(const int n, matrix &unitMat, const ring r=currRing)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
void matrixBlock(const matrix aMat, const matrix bMat, matrix &block)
Creates a new matrix which contains the first argument in the top left corner, and the second in the ...
void henselFactors(const int xIndex, const int yIndex, const poly h, const poly f0, const poly g0, const int d, poly &f, poly &g)
Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a certain degree in x,...
bool pivot(const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring r=currRing)
Finds the best pivot element in some part of a given matrix.
int similar(const number *nn, const int nnLength, const number n, const number tolerance)
Tries to find the number n in the array nn[0..nnLength-1].
bool qrDS(const int n, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
void lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)
LU-decomposition of a given (m x n)-matrix with performing only those divisions that yield zero remai...
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
int luRank(const matrix aMat, const bool isRowEchelon, const ring r=currRing)
Computes the rank of a given (m x n)-matrix.
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.
bool subMatrix(const matrix aMat, const int rowIndex1, const int rowIndex2, const int colIndex1, const int colIndex2, matrix &subMat)
Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.
int pivotScore(number n, const ring r=currRing)
Returns a pivot score for any non-zero matrix entry.
bool upperRightTriangleInverse(const matrix uMat, matrix &iMat, bool diagonalIsOne, const ring r=currRing)
Computes the inverse of a given (n x n)-matrix in upper right triangular form.
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define block
Definition: scanner.cc:666
#define R
Definition: sirandom.c:27