Changeset 12cca3 in git


Ignore:
Timestamp:
May 28, 2010, 2:56:59 PM (14 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
6a649d9e43f6a5e7fd3f1195533f57fb0e5dfccc
Parents:
747a1eef38fad1550aabd46be9267c9003769ac9
Message:
new commands ludecomp(mat), and luinverse(mat [, mat, mat]) (see linearAlgebra.h for docu)

git-svn-id: file:///usr/local/Singular/svn/trunk@12828 2c84dea3-7e68-4137-9b89-c4e89433aadc
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • Singular/iparith.cc

    r747a1e r12cca3  
    44564456  matrix pMat;
    44574457  matrix lMat;
    4458   matrix uMat = mpNew(rr, cc);
    4459 
    4460   /* make initial settings: */
    4461   for (int r = 1; r <= rr; r++)
    4462     for (int c = 1; c <= cc; c++)
    4463       /* 'rMat' is a copy of 'mat': */
    4464       MATELEM(uMat, r, c) = pCopy(mat->m[c - 1 + (r - 1) * cc]);
    4465 
    4466   luDecomp(pMat, lMat, uMat);
     4458  matrix uMat;
     4459
     4460  luDecomp(mat, pMat, lMat, uMat);
    44674461
    44684462  lists ll = (lists)omAllocBin(slists_bin);

  • kernel/LinearAlgebra.cc

    r747a1e r12cca3  
     1/*****************************************************************************\
     2 * Computer Algebra System SINGULAR   
     3\*****************************************************************************/
     4/** @file lineareAlgebra.cc
     5 *
     6 * This file implements basic linear algebra functionality.
     7 *
     8 * For more general information, see the documentation in
     9 * lineareAlgebra.h.
     10 *
     11 * @author Frank Seelisch
     12 *
     13 * @internal @version \$Id$
     14 *
     15 **/
     16/*****************************************************************************/
     17
     18// include header files
    119#include "mod2.h"
    220#include "structs.h"
    … …  
    523#include "numbers.h"
    624#include "matpol.h"
    7 #include "LinearAlgebra.h"
    8 
    9 int pivotStrategy(const number n)
    10 { /* bear in mind that n is guaranteed to be non-zero */
    11   return 0;
    12  
    13   /* in K(x): degree(denominator) - degree(numerator) -> maximal!
    14      in Q:    abs(n) -> maximal
    15      sonst:   0 */
    16 }
    17 
     25#include "linearAlgebra.h"
     26
     27/**
     28 * The returned score is based on the implementation of 'nSize' for
     29 * numbers (, see numbers.h): nSize(n) provides a measure for the
     30 * complexity of n. Thus, less complex pivot elements will be
     31 * prefered, and get therefore a smaller pivot score. Consequently,
     32 * we simply return the value of nSize.
     33 * An exception to this rule are the ground fields R, long R, and
     34 * long C: In these, the result of nSize relates to |n|. And since
     35 * a larger modulus of the pivot element ensures a numerically more
     36 * stable Gauss elimination, we would like to have a smaller score
     37 * for larger values of |n|. Therefore, in R, long R, and long C,
     38 * the negative of nSize will be returned.
     39 **/
     40int pivotScore(number n)
     41{
     42  int s = nSize(n);
     43  if (rField_is_long_C(currRing) ||
     44      rField_is_long_R(currRing) ||
     45      rField_is_R(currRing))
     46    return -s;
     47  else
     48    return s;
     49}
     50
     51/**
     52 * This code computes a score for each non-zero matrix entry in
     53 * aMat[r1..r2, c1..c2]. If all entries are zero, false is returned,
     54 * otherwise true. In the latter case, the minimum of all scores
     55 * is sought, and the row and column indices of the corresponding
     56 * matrix entry are stored in bestR and bestC.
     57 **/
    1858bool pivot(const matrix aMat, const int r1, const int r2, const int c1,
    1959           const int c2, int* bestR, int* bestC)
    2060{
    21   int bestScore = -1; int score; poly matEntry;
     61  int bestScore; int score; bool foundBestScore = false; poly matEntry;
    2262
    2363  for (int c = c1; c <= c2; c++)
    … …  
    2666    {
    2767      matEntry = MATELEM(aMat, r, c);
    28       score = (matEntry == NULL) ? -1 : pivotStrategy(pGetCoeff(matEntry));
    29       if (score > bestScore)
    30       {
    31         bestScore = score;
    32         *bestR = r;
    33         *bestC = c;
    34       }
    35     }
    36   }
    37 
    38   return (bestScore != -1);
    39 }
    40 
    41 void luDecomp(matrix &pMat, matrix &lMat, matrix uMat)
    42 {
    43   int rr = uMat->rows();
    44   int cc = uMat->cols();
     68      if (matEntry != NULL)
     69      {
     70        score = pivotScore(pGetCoeff(matEntry));
     71        if ((!foundBestScore) || (score < bestScore))
     72        {
     73          bestScore = score;
     74          *bestR = r;
     75          *bestC = c;
     76        }
     77        foundBestScore = true;
     78      }
     79    }
     80  }
     81
     82  return foundBestScore;
     83}
     84
     85void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
     86{
     87  int rr = aMat->rows();
     88  int cc = aMat->cols();
    4589  pMat = mpNew(rr, rr);
    4690  lMat = mpNew(rr, rr);
     91  uMat = mpNew(rr, cc);
     92  /* copy aMat into uMat: */
     93  for (int r = 1; r <= rr; r++)
     94    for (int c = 1; c <= cc; c++)
     95      MATELEM(uMat, r, c) = pCopy(aMat->m[c - 1 + (r - 1) * cc]);
     96
    4797  /* we use an int array to store all row permutations;
    4898     note that we only make use of the entries [1..rr] */
    … …  
    73123
    74124      /* swap rows with indices r and bestR in lMat;
    75          we must do this for columns < r */
     125         we must do this only for columns < r */
    76126      for (int c = 1; c < r; c++)
    77127      {
    … …  
    82132
    83133      /* perform next Gauss elimination step, i.e., below the
    84          row with index r; we need to adjust lMat and uMat;
    85          we know that the matrix entry at [r, r] is non-zero: */
     134         row with index r;
     135         we need to adjust lMat and uMat;
     136         we are certain that the matrix entry at [r, r] is non-zero: */
    86137      number pivotElement = pGetCoeff(MATELEM(uMat, r, r));
    87138      poly p; number n;
    … …  
    101152          n = nNeg(n);
    102153          for (int cGauss = r + 1; cGauss <= cc; cGauss++)
    103             MATELEM(uMat, rGauss, cGauss) = pAdd(MATELEM(uMat, rGauss, cGauss),
    104                                                      ppMult_nn(MATELEM(uMat, r,
    105                                                                  cGauss), n));
     154            MATELEM(uMat, rGauss, cGauss)
     155              = pAdd(MATELEM(uMat, rGauss, cGauss),
     156                     ppMult_nn(MATELEM(uMat, r, cGauss), n));
    106157
    107158          nDelete(&n); /* clean up n */
    … …  
    119170}
    120171
     172/**
     173 * This code first computes the LU-decomposition of aMat,
     174 * and then calls the method for inverting a matrix which
     175 * is given by its LU-decomposition.
     176 **/
    121177bool luInverse(const matrix aMat, matrix &iMat)
    122 { /* aMat is guaranteed to be quadratic */
     178{ /* aMat is guaranteed to be an (n x n)-matrix */
    123179  int d = aMat->rows();
    124180
    125   /* preparing to call luInverseFromLUDecomp */
    126   matrix pMat = mpNew(d, d);
    127   matrix lMat = mpNew(d, d);
    128   matrix uMat = mpNew(d, d);
    129   /* copy aMat to uMat: */
    130   for (int r = 1; r <= d; r++)
    131     for (int c = 1; c <= d; c++)
    132       MATELEM(uMat, r, c) = pCopy(aMat->m[c - 1 + (r - 1) * d]);
    133 
    134   luDecomp(pMat, lMat, uMat);
     181  matrix pMat;
     182  matrix lMat;
     183  matrix uMat;
     184  luDecomp(aMat, pMat, lMat, uMat);
    135185  bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat);
    136186
    … …  
    238288}
    239289
    240 bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat)
     290/**
     291 * This code computes the inverse by inverting lMat and uMat, and
     292 * then performing two matrix multiplications.
     293 **/
     294bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
     295                           const matrix uMat, matrix &iMat)
    241296{ /* uMat is guaranteed to be quadratic */
    242297  int d = uMat->rows();
    243298
    244   matrix lMinus1Mat; /* for storing the inverse of lMat */
    245   matrix uMinus1Mat; /* for storing the inverse of uMat, if existing */
    246 
    247   bool result = upperRightTriangleInverse(uMat, uMinus1Mat, false);
     299  matrix lMatInverse; /* for storing the inverse of lMat;
     300                         this inversion will always work                */
     301  matrix uMatInverse; /* for storing the inverse of uMat, if invertible */
     302
     303  bool result = upperRightTriangleInverse(uMat, uMatInverse, false);
    248304  if (result)
    249305  {
    250     lowerLeftTriangleInverse(lMat, lMinus1Mat, true);
    251     iMat = mpMult(mpMult(uMinus1Mat, lMinus1Mat), pMat);
     306    /* next will always work, since lMat is known to have all diagonal
     307       entries equal to 1 */
     308    lowerLeftTriangleInverse(lMat, lMatInverse, true);
     309    iMat = mpMult(mpMult(uMatInverse, lMatInverse), pMat);
    252310   
    253311    /* clean-up */
    254     idDelete((ideal*)&lMinus1Mat);
    255     idDelete((ideal*)&uMinus1Mat);
     312    idDelete((ideal*)&lMatInverse);
     313    idDelete((ideal*)&uMatInverse);
    256314  }
    257315

  • kernel/LinearAlgebra.h

    r747a1e r12cca3  
     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 * @internal @version \$Id$
     20 *
     21 **/
     22/*****************************************************************************/
     23
    124#ifndef LINEAR_ALGEBRA_H
    225#define LINEAR_ALGEBRA_H
    326
    4 /****************************************
    5 *  Computer Algebra System SINGULAR     *
    6 ****************************************/
    7 /* $Id: LineareAlgebra.h$ */
    8 /*
    9  * ABSTRACT: linear algebra over any SINGULAR-supported field
    10  *
    11  * Use 'K' in all method documentations to denote the field.
    12  * All matrices herein are assumed to have 'number' entries,
    13  * representing elements of K.
    14  * Use 'r' and 'c' as canonical row and column indices to ease
    15  * readability of the code.
    16  */
    17 
     27// include basic SINGULAR structures
    1828#include <structs.h>
    1929
    2030/**
    21   * LU-decomposition of a given (m x n)-matrix.
    22   *
    23   * Given an (m x n) matrix A, the method computes its LU-decomposition,
    24   * that is, it computes matrices P, L, and U such that
    25   * - A = P * L * U,
    26   * - P is an (m x m) permutation matrix, i.e., its row/columns are the
    27   *   standard basis vectors of K^m
    28   * - L is an (m x m) matrix in lower triangular form,
    29   * - U is an (m x n) matrix in upper row echelon form
    30   * From these conditions, it easily follows that also
    31   * P * A = L * U, as P is self-inverse.
    32   * The method assumes that U contains the matrix A for which the
    33   * decomposition is being sought.
    34   * The method will modify all agument matrices so that they will
    35   * finally contain the matrices P, L, and U with the above properties.
    36   *
    37   * @param pMat afterwards the row permutation matrix P
    38   * @param lMat afterwards the lower triangular matrix L
    39   * @param uMat afterwards the upper row echelon matrix U, initially A
    40   */
    41 void luDecomp(matrix &pMat, matrix &lMat, matrix uMat);
    42 
    43 /**
    44   * Finds the best pivot element in some part of a given matrix.
    45   *
    46   * Given any matrix A with valid row indices r1..r2 and valid column
    47   * indices c1..c2, this method finds the best pivot element for
    48   * continuing Gauss elimination in A. "Best" here means best according
    49   * to the implemented pivotStrategy; see there.
    50   * In case all elements in A[r1..r2, c1..c2] are zero, the method returns
    51   * false, otherwise true.
    52   *
    53   * @param aMat any matrix
    54   * @param r1 the lower row index, determines the relevant part of aMat
    55   * @param r2 the upper row index, determines the relevant part of aMat
    56   * @param c1 the lower column index, determines the relevant part of aMat
    57   * @param c2 the upper column index, determines the relevant part of aMat
    58   * @param bestR afterwards address of the row index of the best pivot element
    59   * @param bestC afterwards address of the col index of the best pivot element
    60   * @return false if all relevant matrix entries are zero, false otherwise
    61   */
    62 bool pivot(const matrix aMat, const int r1, const int r2, const int c1,
    63            const int c2, int* bestR, int* bestC);
    64 
    65 /**
    66   * Computes an estimate for how well a given non-zero number can serve as
    67   * pivot element for the next Gauss elimination step.
    68   *
    69   * This method expects a non-zero number. The following rules must be obeyed:
    70   * - return a non-negative int
    71   * - the higher the value, the better will n serve as pivot element
    72   *
    73   * @param n a number representing a non-zero element of the ground field K
    74   * @return estimate for how well n will serve as pivot element
    75   */
    76 int pivotStrategy(const number n);
    77 
    78 /**
    79   * Computes the inverse of a quadratic matrix.
    80   *
    81   * This method expects a quadratic matrix, that is, it must have as many rows
    82   * as columns. Inversion of the first argument is attempted via the
    83   * LU-decomposition. There are two cases:
    84   * 1) The matrix is not invertible. Then the method returns false, and the
    85   *    content of iMat is useless.
    86   * 2) The matrix is invertible. Then the method returns true, and the
    87   *    content of iMat is the inverse of aMat.
    88   *
    89   * @param aMat the quadratic matrix to be inverted
    90   * @param afterwards iMat the inverse of aMat in case that aMat is invertible
    91   * @return true iff aMat is invertible, false otherwise
    92   */
    93 bool luInverse(const matrix aMat, matrix &iMat);
    94 
    95 /**
    96   * Computes the inverse of a quadratic matrix in upper right row echelon form.
    97   *
    98   * This method expects a quadratic matrix, that is, it must have as many rows
    99   * as columns. Moreover, uMat[i,j] = 0, at least for all i > j.
    100   * If the argument diagonalIsOne is true, then we know additionally, that
    101   * uMat[i,i] = 1, for all i. In this case uMat is invertible. If diagonalIsOne
    102   * is false, we do not know anything about the diagonal entries. (Note that
    103   * they may still all be 1.)
    104   * In general, there are two cases:
    105   * 1) The matrix is not invertible. Then the method returns false, and the
    106   *    content of iMat is useless.
    107   * 2) The matrix is invertible. Then the method returns true, and the
    108   *    content of iMat is the inverse of aMat.
    109   *
    110   * @param uMat a quadratic matrix in upper right row echelon form
    111   * @param iMat afterwards the inverse of uMat in case that uMat is invertible
    112   * @param diagonalIsOne if true all diagonal entries of uMat are 1
    113   * @return true iff uMat is invertible, false otherwise
    114   */
    115 bool upperRightTriangleInverse(const matrix uMat, matrix &iMat,
    116                                bool diagonalIsOne);
    117 
    118 /**
    119   * Computes the inverse of a quadratic matrix in lower left row echelon form.
    120   *
    121   * This method expects a quadratic matrix, that is, it must have as many rows
    122   * as columns. Moreover, lMat[i,j] = 0, at least for all j > i.
    123   * If the argument diagonalIsOne is true, then we know additionally, that
    124   * lMat[i,i] = 1, for all i. In this case lMat is invertible. If diagonalIsOne
    125   * is false, we do not know anything about the diagonal entries. (Note that
    126   * they may still all be 1.)
    127   * In general, there are two cases:
    128   * 1) The matrix is not invertible. Then the method returns false, and the
    129   *    content of iMat is useless.
    130   * 2) The matrix is invertible. Then the method returns true, and the
    131   *    content of iMat is the inverse of aMat.
    132   *
    133   * @param lMat a quadratic matrix in lower left row echelon form
    134   * @param iMat afterwards the inverse of lMat in case that lMat is invertible
    135   * @param diagonalIsOne if true all diagonal entries of lMat are 1
    136   * @return true iff lMat is invertible, false otherwise
    137   */
    138 bool lowerLeftTriangleInverse(const matrix uMat, matrix &iMat,
    139                               bool diagonalIsOne);
    140 
    141 /**
    142   * Computes the inverse of a quadratic matrix which is given by its LU-
    143   * decomposition.
    144   *
    145   * With A denoting the matrix to be inverted, the method expects the
    146   * LU-decomposition of A, that is, pMat * A = lMat * uMat, where
    147   * the argument matrices have the appropriate proteries.
    148   * Furthermore, uMat is expected to be quadratic. Then A^(-1) is computed
    149   * as A^(-1) = uMat^(-1) * lMat^(-1) * pMat, since pMat is self-inverse.
    150   * This will work if and only if uMat is invertible, as lMat and pMat are
    151   * in any case invertible.
    152   *
    153   * There are two cases:
    154   * 1) A is not invertible. Then the method returns false, and the
    155   *    content of iMat is useless.
    156   * 2) A is invertible. Then the method returns true, and the content
    157   *    of iMat is the inverse of A.
    158   *
    159   * @param pMat the permutation matrix of the LU-decomposition of A
    160   * @param lMat the lower triangle matrix of the LU-decomposition of A
    161   * @param uMat the upper row echelon matrix of the LU-decomposition of A
    162   * @param iMat afterwards the inverse of A in case that A is invertible
    163   * @return true iff A is invertible, false otherwise
    164   */
    165 bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
    166                            const matrix uMat, matrix &iMat);
     31 * LU-decomposition of a given (m x n)-matrix.
     32 *
     33 * Given an (m x n) matrix A, the method computes its LU-decomposition,
     34 * that is, it computes matrices P, L, and U such that<br>
     35 * - P * A = L * U,<br>
     36 * - P is an (m x m) permutation matrix, i.e., its row/columns form the
     37 *   standard basis of K^m,<br>
     38 * - L is an (m x m) matrix in lower triangular form with all diagonal
     39 *   entries equal to 1,<br>
     40 * - U is an (m x n) matrix in upper row echelon form.<br>
     41 * From these conditions, it easily follows that also A = P * L * U,
     42 * since P is self-inverse.<br>
     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             );
     53
     54/**
     55 * Returns a pivot score for any non-zero matrix entry.
     56 *
     57 * The smaller the score the better will n serve as a pivot element
     58 * in subsequent Gauss elimination steps.
     59 *
     60 * @return the pivot score of n
     61 **/
     62int pivotScore(
     63               number n  /**< [in] a non-zero matrix entry */
     64              );
     65
     66/**
     67 * Finds the best pivot element in some part of a given matrix.
     68 *
     69 * Given any matrix A with valid row indices r1..r2 and valid column
     70 * indices c1..c2, this method finds the best pivot element for
     71 * subsequent Gauss elimination steps in A[r1..r2, c1..c2]. "Best"
     72 * here means best with respect to the implementation of the method
     73 * 'pivotScore(number n)'.<br>
     74 * In the case that all elements in A[r1..r2, c1..c2] are zero, the
     75 * method returns false, otherwise true.
     76 *
     77 * @return false if all relevant matrix entries are zero, true otherwise
     78 * @sa pivotScore(number n)
     79 **/
     80bool pivot(
     81           const matrix aMat, /**< [in]  any matrix with number entries */
     82           const int r1,      /**< [in]  lower row index                */
     83           const int r2,      /**< [in]  upper row index                */
     84           const int c1,      /**< [in]  lower column index             */
     85           const int c2,      /**< [in]  upper column index             */
     86           int* bestR,        /**< [out] address of row index of best
     87                                         pivot element                  */
     88           int* bestC         /**< [out] address of column index of
     89                                         best pivot element             */
     90          );
     91
     92/**
     93 * Computes the inverse of a given (n x n)-matrix.
     94 *
     95 * This method expects an (n x n)-matrix, that is, it must have as many
     96 * rows as columns. Inversion of the first argument is attempted via the
     97 * LU-decomposition. There are two cases:<br>
     98 * 1) The matrix is not invertible. Then the method returns false, and
     99 *    &iMat remains unchanged.<br>
     100 * 2) The matrix is invertible. Then the method returns true, and the
     101 *    content of iMat is the inverse of aMat.
     102 *
     103 * @return true iff aMat is invertible, false otherwise
     104 * @sa luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
     105 *                           const matrix uMat, matrix &iMat)
     106 **/
     107bool luInverse(
     108               const matrix aMat, /**< [in]  matrix to be inverted */
     109               matrix &iMat       /**< [out] inverse of aMat if
     110                                             invertible            */
     111              );
     112
     113/**
     114 * Computes the inverse of a given (n x n)-matrix in upper right
     115 * triangular form.
     116 *
     117 * This method expects a quadratic matrix, that is, it must have as
     118 * many rows as columns. Moreover, uMat[i, j] = 0, at least for all
     119 * i > j, that is, u is in upper right triangular form.<br>
     120 * If the argument diagonalIsOne is true, then we know additionally,
     121 * that uMat[i, i] = 1, for all i. In this case uMat is invertible.
     122 * Contrariwise, if diagonalIsOne is false, we do not know anything
     123 * about the diagonal entries. (Note that they may still all be
     124 * 1.)<br>
     125 * In general, there are two cases:<br>
     126 * 1) The matrix is not invertible. Then the method returns false,
     127 *    and &iMat remains unchanged.<br>
     128 * 2) The matrix is invertible. Then the method returns true, and
     129 *    the content of iMat is the inverse of uMat.
     130 *
     131 * @return true iff uMat is invertible, false otherwise
     132 **/
     133bool upperRightTriangleInverse(
     134       const matrix uMat, /**< [in]  (n x n)-matrix in upper right
     135                                     triangular form               */
     136       matrix &iMat,      /**< [out] inverse of uMat if invertible */
     137       bool diagonalIsOne /**< [in]  if true, then all diagonal
     138                                     entries of uMat are 1         */
     139                              );
     140
     141/**
     142 * Computes the inverse of a given (n x n)-matrix in lower left
     143 * triangular form.
     144 *
     145 * This method expects an (n x n)-matrix, that is, it must have as
     146 * many rows as columns. Moreover, lMat[i,j] = 0, at least for all
     147 * j > i, that ism lMat is in lower left triangular form.<br>
     148 * If the argument diagonalIsOne is true, then we know additionally,
     149 * that lMat[i, i] = 1, for all i. In this case lMat is invertible.
     150 * Contrariwise, if diagonalIsOne is false, we do not know anything
     151 * about the diagonal entries. (Note that they may still all be
     152 * 1.)<br>
     153 * In general, there are two cases:<br>
     154 * 1) The matrix is not invertible. Then the method returns false,
     155 *    and &iMat remains unchanged.<br>
     156 * 2) The matrix is invertible. Then the method returns true, and
     157 *    the content of iMat is the inverse of lMat.
     158 *
     159 * @return true iff lMat is invertible, false otherwise
     160 **/
     161bool lowerLeftTriangleInverse(
     162       const matrix lMat, /**< [in]  (n x n)-matrix in lower left
     163                                     triangular form               */
     164       matrix &iMat,      /**< [out] inverse of lMat if invertible */
     165       bool diagonalIsOne /**< [in]  if true, then all diagonal
     166                                     entries of lMat are 1         */
     167                              );
     168
     169/**
     170 * Computes the inverse of an (n x n)-matrix which is given by its LU-
     171 * decomposition.
     172 *
     173 * With A denoting the matrix to be inverted, the method expects the
     174 * LU-decomposition of A, that is, pMat * A = lMat * uMat, where
     175 * the argument matrices have the appropriate proteries; see method
     176 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
     177 * matrix &uMat)'.<br>
     178 * Furthermore, uMat is expected to be an (n x n)-matrix. Then A^(-1)
     179 * is computed according to A^(-1) = uMat^(-1) * lMat^(-1) * pMat,
     180 * since pMat is self-inverse. This will work if and only if uMat is
     181 * invertible, because lMat and pMat are in any case invertible.<br>
     182 * In general, there are two cases:<br>
     183 * 1) uMat and hence A is not invertible. Then the method returns
     184 *    false, and &iMat remains unchanged.<br>
     185 * 2) uMat and hence A is invertible. Then the method returns true,
     186 *    and the content of iMat is the inverse of A.
     187 *
     188 * @return true if A is invertible, false otherwise
     189 * @sa luInverse(const matrix aMat, matrix &iMat)
     190 **/
     191bool luInverseFromLUDecomp(
     192       const matrix pMat, /**< [in]  permutation matrix of an LU-
     193                                     decomposition                */
     194       const matrix lMat, /**< [in]  lower left matrix of an LU-
     195                                     decomposition                */
     196       const matrix uMat, /**< [in]  upper right matrix of an LU-
     197                                     decomposition                */
     198       matrix &iMat       /**< [out] inverse of A if invertible   */
     199                          );
    167200
    168201#endif

  • kernel/Makefile.in

    r747a1e r12cca3  
    115115    pShallowCopyDelete.cc fast_mult.cc digitech.cc \
    116116    tgb.cc tgbgauss.cc ringgb.cc f5data.cc f5lists.cc f5gb.cc \
    117     f5c.cc F5cLists.cc ratgring.cc shiftgb.cc gfan.cc LinearAlgebra.cc
     117    f5c.cc F5cLists.cc ratgring.cc shiftgb.cc gfan.cc linearAlgebra.cc
    118118
    119119# normal C source files

  • kernel/gnumpc.cc

    r747a1e r12cca3  
    120120}
    121121
     122int ngcSize(number n)
     123{
     124  int r = (int)((gmp_complex*)n)->real();
     125  if (r < 0) r = -r;
     126  int i = (int)((gmp_complex*)n)->imag();
     127  if (i < 0) i = -i;
     128  int oneNorm = r + i;
     129  /* basically return the 1-norm of n;
     130     only if this happens to be zero although n != 0,
     131     return 1;
     132     (this code ensures that zero has the size zero) */
     133  if ((oneNorm == 0.0) & (ngcIsZero(n) == FALSE)) oneNorm = 1;
     134  return oneNorm;
     135}
     136
    122137/*2
    123138* delete a

  • kernel/gnumpc.h

    r747a1e r12cca3  
    2525number   ngcMult(number a, number b);
    2626number   ngcDiv(number a, number b);
     27int      ngcSize(number n);
    2728void     ngcPower(number x, int exp, number *lu);
    2829number   ngcCopy(number a);

  • kernel/gnumpfl.cc

    r747a1e r12cca3  
    9797  else
    9898    return (int)(d+0.5);
     99}
     100
     101int ngfSize(number n)
     102{
     103  int i = ngfInt(n, currRing);
     104  /* basically return the largest integer in n;
     105     only if this happens to be zero although n != 0,
     106     return 1;
     107     (this code ensures that zero has the size zero) */
     108  if ((i == 0) && (ngfIsZero(n) == FALSE)) i = 1;
     109  return i;
    99110}
    100111

  • kernel/gnumpfl.h

    r747a1e r12cca3  
    2323number   ngfSub(number la, number li);
    2424number   ngfMult(number a, number b);
     25int      ngfSize(number n);
    2526number   ngfDiv(number a, number b);
    2627void     ngfPower(number x, int exp, number *lu);

  • kernel/numbers.cc

    r747a1e r12cca3  
    607607    n->cfSetMap=nrSetMap;
    608608    /* nName= ndName; */
    609     /*nSize  = ndSize;*/
     609    n->nSize = nrSize;
    610610#ifdef LDEBUG
    611611    n->nDBTest=ndDBTest; // not yet implemented: nrDBTest;
    … …  
    637637    n->cfSetMap=ngfSetMap;
    638638    n->nName= ndName;
    639     n->nSize  = ndSize;
     639    n->nSize  = ngfSize;
    640640#ifdef LDEBUG
    641641    n->nDBTest=ndDBTest; // not yet implemented: ngfDBTest
    … …  
    670670    n->nRePart=ngcRePart;
    671671    n->nImPart=ngcImPart;
    672     /*nSize  = ndSize;*/
     672    n->nSize = ngcSize;
    673673#ifdef LDEBUG
    674674    n->nDBTest=ndDBTest; // not yet implemented: ngcDBTest

  • kernel/numbers.h

    r747a1e r12cca3  
    5151extern number  (*nPar)(int i);
    5252extern int     (*nParDeg)(number n);
     53/* size: a measure for the complexity of the represented number n;
     54         zero should have size zero; larger size means more complex */
    5355extern int     (*nSize)(number n);
    5456extern int     (*n_Int)(number &n, const ring r);

  • kernel/shortfl.cc

    r747a1e r12cca3  
    6565  else
    6666    i = 0;
     67  return i;
     68}
     69
     70int nrSize(number n)
     71{
     72  float f = nf(n).F();
     73  int i = (int)f;
     74  /* basically return the largest integer in n;
     75     only if this happens to be zero although n != 0,
     76     return 1;
     77     (this code ensures that zero has the size zero) */
     78  if ((f != 0.0) & (i == 0)) i = 1;
    6779  return i;
    6880}

  • kernel/shortfl.h

    r747a1e r12cca3  
    2323number  nrNeg         (number c);
    2424number  nrInvers      (number c);
     25int     nrSize        (number n);
    2526BOOLEAN nrGreater     (number a, number b);
    2627BOOLEAN nrEqual       (number a, number b);
Note: See TracChangeset for help on using the changeset viewer.