My Project  debian-1:4.1.1-p2+ds-4build3
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 */
subMatrix
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.
Definition: linearAlgebra.cc:733
pivot
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.
Definition: linearAlgebra.cc:68
swapRows
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
Definition: linearAlgebra.cc:781
charPoly
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
Definition: linearAlgebra.cc:748
ip_smatrix
Definition: matpol.h:15
f
FILE * f
Definition: checklibs.c:9
pivotScore
int pivotScore(number n, const ring r=currRing)
Returns a pivot score for any non-zero matrix entry.
Definition: linearAlgebra.cc:50
luRank
int luRank(const matrix aMat, const bool isRowEchelon, const ring r=currRing)
Computes the rank of a given (m x n)-matrix.
Definition: linearAlgebra.cc:230
quadraticSolve
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
Definition: linearAlgebra.cc:626
lduDecomp
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...
Definition: linearAlgebra.cc:1343
g
g
Definition: cfModGcd.cc:4031
realSqrt
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.
Definition: linearAlgebra.cc:601
luDecomp
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring r=currRing)
LU-decomposition of a given (m x n)-matrix.
Definition: linearAlgebra.cc:103
currRing
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
matrixBlock
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 ...
Definition: linearAlgebra.cc:805
henselFactors
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,...
Definition: linearAlgebra.cc:1219
block
#define block
Definition: scanner.cc:665
structs.h
h
static Poly * h
Definition: janet.cc:972
upperRightTriangleInverse
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.
Definition: linearAlgebra.cc:251
luSolveViaLUDecomp
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...
Definition: linearAlgebra.cc:377
luInverse
bool luInverse(const matrix aMat, matrix &iMat, const ring r=currRing)
Computes the inverse of a given (n x n)-matrix.
Definition: linearAlgebra.cc:200
luInverseFromLUDecomp
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.
Definition: linearAlgebra.cc:352
unitMatrix
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.
Definition: linearAlgebra.cc:95
lowerLeftTriangleInverse
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.
Definition: linearAlgebra.cc:300
swapColumns
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
Definition: linearAlgebra.cc:793
l
int l
Definition: cfEzgcd.cc:93
hessenberg
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring r)
Computes the Hessenberg form of a given square matrix.
Definition: linearAlgebra.cc:909
R
#define R
Definition: sirandom.c:26
p
int p
Definition: cfModGcd.cc:4019
H
CanonicalForm H
Definition: facAbsFact.cc:64
qrDS
bool qrDS(const int n, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
Definition: linearAlgebra.cc:1090
similar
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].
Definition: linearAlgebra.cc:1188
luSolveViaLDUDecomp
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...
Definition: linearAlgebra.cc:1461