AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MA28_CppDecl.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 //
00042 // C declarations for MA28 functions.  These declarations should not have to change
00043 // for different platforms.  As long as the fortran object code uses capitalized
00044 // names for its identifers then the declarations in Teuchos_F77_wrappers.h should be
00045 // sufficent for portability.
00046 
00047 #ifndef MA28_CPPDECL_H
00048 #define MA28_CPPDECL_H
00049 
00050 #include "Teuchos_F77_wrappers.h"
00051 
00052 namespace MA28_CppDecl {
00053 
00054 // Declarations that will link to the fortran object file.
00055 // These may change for different platforms
00056 
00057 using FortranTypes::f_int;      // INTEGER
00058 using FortranTypes::f_real;     // REAL
00059 using FortranTypes::f_dbl_prec;   // DOUBLE PRECISION
00060 using FortranTypes::f_logical;    // LOGICAL
00061 
00062 namespace Fortran {
00063 extern "C" {
00064 
00065 // analyze and factorize a matrix
00066 FORTRAN_FUNC_DECL_UL(void,MA28AD,ma28ad) (const f_int& n, const f_int& nz, f_dbl_prec a[], const f_int& licn
00067   , f_int irn[], const f_int& lirn, f_int icn[], const f_dbl_prec& u, f_int ikeep[], f_int iw[]
00068   , f_dbl_prec w[], f_int* iflag);
00069   
00070 // factor using previous analyze
00071 FORTRAN_FUNC_DECL_UL(void,MA28BD,ma28bd) (const f_int& n, const f_int& nz, f_dbl_prec a[], const f_int& licn
00072   , const f_int ivect[], const f_int jvect[], const f_int icn[], const f_int ikeep[], f_int iw[]
00073   , f_dbl_prec w[], f_int* iflag);
00074 
00075 // solve for rhs using internally stored factorized matrix
00076 FORTRAN_FUNC_DECL_UL(void,MA28CD,ma28cd) (const f_int& n, const f_dbl_prec a[], const f_int& licn, const f_int icn[]
00077   , const f_int ikeep[], f_dbl_prec rhs[], f_dbl_prec w[], const f_int& mtype);
00078 
00079 // /////////////////////////////////////////////////////////////////////////////////////////
00080 // Declare structs that represent the MA28 common blocks.  
00081 // These are the common block variables that are ment to be accessed by the user
00082 // Some are used to set the options of MA28 and others return information
00083 // about the attempts to solve the system.
00084 // I want to provide the access functions that allow all of those common block
00085 // variables that are ment to be accessed by the user to be accessable.
00086 // For each of the common data items there will be a get operation that 
00087 // returns the variable value.  For those items that are ment to be
00088 // set by the user there will also be set operations.
00089 
00090 //  COMMON /MA28ED/ LP, MP, LBLOCK, GROW
00091 //  INTEGER LP, MP
00092 //  LOGICAL LBLOCK, GROW
00093 struct MA28ED_struct {
00094   f_int   lp;
00095   f_int   mp;
00096   f_logical lblock;
00097   f_logical grow;
00098 };
00099 extern MA28ED_struct FORTRAN_NAME_UL(MA28ED,ma28ed); // link to fortan common block
00100 
00101 //  COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
00102 // * IRANK, ABORT1, ABORT2
00103 //  INTEGER IRNCP, ICNCP, MINIRN, MINICN, IRANK
00104 //  LOGICAL ABORT1, ABORT2
00105 //  REAL EPS, RMIN, RESID
00106 struct MA28FD_struct {
00107   f_dbl_prec  eps;
00108   f_dbl_prec  rmin;
00109   f_dbl_prec  resid;
00110   f_int   irncp;
00111   f_int   icncp;
00112   f_int   minirn;
00113   f_int   minicn;
00114   f_int   irank;
00115   f_logical abort1;
00116   f_logical abort2;
00117 };
00118 extern MA28FD_struct FORTRAN_NAME_UL(MA28FD,ma28fd); // link to fortan common block
00119 
00120 
00121 //  COMMON /MA28GD/ IDISP
00122 //  INTEGER IDISP
00123 struct MA28GD_struct {
00124   f_int   idisp[2];
00125 };
00126 extern MA28GD_struct FORTRAN_NAME_UL(MA28GD,ma28gd); // link to fortan common block
00127 
00128 //  COMMON /MA28HD/ TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE,
00129 // * NDROP, MAXIT, NOITER, NSRCH, ISTART, LBIG
00130 //  INTEGER NDROP, MAXIT, NOITER, NSRCH, ISTART
00131 //  LOGICAL LBIG
00132 //  REAL TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE
00133 struct MA28HD_struct {
00134   f_dbl_prec  tol;
00135   f_dbl_prec  themax;
00136   f_dbl_prec  big;
00137   f_dbl_prec  dxmax;
00138   f_dbl_prec  errmax;
00139   f_dbl_prec  dres;
00140   f_dbl_prec  cgce;
00141   f_int   ndrop;
00142   f_int   maxit;
00143   f_int   noiter;
00144   f_int   nsrch;
00145   f_int   istart;
00146   f_logical lbig;
00147 };
00148 extern MA28HD_struct FORTRAN_NAME_UL(MA28HD,ma28hd); // link to fortan common block
00149 
00150 //  COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
00151 //  INTEGER LP
00152 //  LOGICAL ABORT1, ABORT2, ABORT3
00153 struct MA30ED_struct {
00154   f_int   lp;
00155   f_logical abort1;
00156   f_logical abort2;
00157   f_logical abort3;
00158 };
00159 extern MA30ED_struct FORTRAN_NAME_UL(MA30ED,ma30ed); // link to fortan common block
00160 
00161 //  COMMON /MA30FD/ IRNCP, ICNCP, IRANK, IRN, ICN
00162 //  INTEGER IRNCP, ICNCP, IRANK, IRN, ICN
00163 struct MA30FD_struct {
00164   f_int   irncp;
00165   f_int   icncp;
00166   f_int   irank;
00167   f_int   minirn;
00168   f_int   minicn;
00169 };
00170 extern MA30FD_struct FORTRAN_NAME_UL(MA30FD,ma30fd); // link to fortan common block
00171 
00172 //  COMMON /MA30GD/ EPS, RMIN
00173 //  DOUBLE PRECISION EPS, RMIN
00174 struct MA30GD_struct {
00175   f_dbl_prec  eps;
00176   f_dbl_prec  rmin;
00177 };
00178 extern MA30GD_struct FORTRAN_NAME_UL(MA30GD,ma30gd); // link to fortan common block
00179 
00180 //  COMMON /MA30HD/ RESID
00181 //  DOUBLE PRECISION RESID
00182 struct MA30HD_struct {
00183   f_dbl_prec  resid;
00184 };
00185 extern MA30HD_struct FORTRAN_NAME_UL(MA30HD,ma30hd); // link to fortan common block
00186 
00187 //  COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
00188 //  INTEGER NDROP, NSRCH
00189 //  LOGICAL LBIG
00190 //  DOUBLE PRECISION TOL, BIG
00191 struct MA30ID_struct {
00192   f_dbl_prec  tol;
00193   f_dbl_prec  big;
00194   f_int   ndrop;
00195   f_int   nsrch;
00196   f_logical lbig;
00197 };
00198 extern MA30ID_struct FORTRAN_NAME_UL(MA30ID,ma30id); // link to fortan common block
00199 
00200 //  COMMON /MC23BD/ LP,NUMNZ,NUM,LARGE,ABORT
00201 //  INTEGER LP, NUMNZ, NUM, LARGE
00202 //  LOGICAL ABORT
00203 struct MC23BD_struct {
00204   f_int   lp;
00205   f_int   numnz;
00206   f_int   num;
00207   f_int   large;
00208   f_logical abort;
00209 };
00210 extern MC23BD_struct FORTRAN_NAME_UL(MC23BD,mc23bd); // link to fortan common block
00211 
00212 
00213 } // end extern "C"
00214 } // end namespace Fortran
00215 
00216 
00217 /* * @name {\bf MA28 C++ Declarations}.
00218   *
00219   * These the C++ declarations for MA28 functions and common block data.
00220   * These declarations will not change for different platforms.
00221   * All of these functions are in the C++ namespace #MA28_CppDecl#.
00222   *
00223   * These functions perform the three phases that are normally associated
00224   * with solving sparse systems of linear equations; analyze (and factorize)
00225   * , factorize, and solve.  The MA28 interface uses a coordinate format
00226   * (aij, i, j) for the sparse matrix.
00227   * 
00228   * There are three interface routienes that perform these steps:
00229   * \begin{description}
00230   * \item[#ma28ad#] Analyzes and factorizes a sparse matrix stored in #a#, #irn#
00231   *   , #icn# and returns the factorized matrix data structures in #a#, #icn#
00232   *   , and #ikeep#.
00233   * \item[#ma28bd#] Factorizes a matrix with the same sparsity structure that was
00234   *   previously factorized by #ma28ad#.  Information about the row and column
00235   *   permutations, fill-in elements etc. from the previous analyze and factorization
00236   *   is passed in the arguments #icn#, and #ikeep#.  The matrix to be
00237   *   factorized is passed in #a#, #ivect#, and #jvect# and the non-zero
00238   *   elements of the factorization are returned in #a#.
00239   * \item[#ma28cd#] Solves for a dense right hand side (rhs) given a matrix factorized
00240   *   by #ma28ad# or #ma28bd#.  The rhs is passed in #rhs# and the solution
00241   *   is returned in #rhs#.  The factorized matrix is passed in by #a#, #icn#
00242   *   and #ikeep#.  The transposed or the non-transposed system can be solved
00243   *   for by passing in #mtype != 1# and #mtype == 1# respectively.
00244   */
00245 
00246 // @{
00247 //    begin MA28 C++ Declarations
00248 
00249 /* * @name {\bf MA28 / MA30 Common Block Access}.
00250   *
00251   * These are references to structures that allow C++ users to set and retrive
00252   * values of the MA28xD and MA30xD common blocks.   Some of the common block
00253   * items listed below for MA28xD are also present in MA30xD.  The control
00254   * parameters (abort1, eps, etc.) for MA28xD are transfered to the equivalent
00255   * common block variables in the #ma28ad# function but not in any of the other
00256   * functions.
00257   *
00258   * The internal states for MA28, MA30, and MC23 are determined by the 
00259   * values in these common block variables as there are no #SAVE# variables
00260   * in any of the functions.  So to use MA28 with more than
00261   * one sparse matrix at a time you just have to keep copies of these
00262   * common block variable for each system and then set them when every
00263   * you want to work with that system agian.  This is very good news.
00264   *
00265   * These common block variables are:
00266   * \begin{description}
00267   * \item[lp, mp]
00268   *   Integer: Used by the subroutine as the unit numbers for its warning
00269   *   and diagnostic messages. Default value for both is 6 (for line
00270   *   printer output). the user can either reset them to a different
00271   *   stream number or suppress the output by setting them to zero.
00272   *   While #lp# directs the output of error diagnostics from the
00273   *   principal subroutines and internally called subroutines, #mp#
00274   *   controls only the output of a message which warns the user that he
00275   *   has input two or more non-zeros a(i), . . ,a(k) with the same row
00276   *   and column indices.  The action taken in this case is to proceed
00277   *   using a numerical value of a(i)+...+a(k). in the absence of other
00278   *   errors, #iflag# will equal -14 on exit.
00279   * \item[lblock]
00280   *   Logical: Controls an option of first
00281   *   preordering the matrix to block lower triangular form (using
00282   *   harwell subroutine mc23a). The preordering is performed if #lblock#
00283   *   is equal to its default value of #true# if #lblock# is set to
00284   *   #false# , the option is not invoked and the space allocated to
00285   *   #ikeep# can be reduced to 4*n+1.
00286   * \item[grow]
00287   *    Logical: If it is left at its default value of
00288   *   #true# , then on return from ma28a/ad or ma28b/bd, w(1) will give
00289   *   an estimate (an upper bound) of the increase in size of elements
00290   *   encountered during the decomposition. If the matrix is well
00291   *   scaled, then a high value for w(1), relative to the largest entry
00292   *   in the input matrix, indicates that the LU decomposition may be
00293   *   inaccurate and the user should be wary of his results and perhaps
00294   *   increase u for subsequent runs.  We would like to emphasise that
00295   *   this value only relates to the accuracy of our LU decomposition
00296   *   and gives no indication as to the singularity of the matrix or the
00297   *   accuracy of the solution.  This upper bound can be a significant
00298   *   overestimate particularly if the matrix is badly scaled. If an
00299   *   accurate value for the growth is required, #lbig# (q.v.) should be
00300   *   set to #true#
00301   * \item[eps, rmin]
00302   *   Double Precision:  If on entry to ma28b/bd, #eps# is less
00303   *   than one, then #rmin# will give the smallest ratio of the pivot to
00304   *   the largest element in the corresponding row of the upper
00305   *   triangular factor thus monitoring the stability of successive
00306   *   factorizations. if rmin becomes very large and w(1) from
00307   *   ma28b/bd is also very large, it may be advisable to perform a
00308   *    new decomposition using ma28a/ad.
00309   * \item[resid]
00310   *   Double Precision:  On exit from ma28c/cd gives the value
00311   *   of the maximum residual over all the equations unsatisfied because
00312   *   of dependency (zero pivots).
00313   * \item[irncp,icncp]
00314   *   Integer:  Monitors the adequacy of "elbow
00315   *   room" in #irn# and #a#/#icn# respectively. If either is quite large (say
00316   *   greater than n/10), it will probably pay to increase the size of
00317   *   the corresponding array for subsequent runs. if either is very low
00318   *   or zero then one can perhaps save storage by reducing the size of
00319   *   the corresponding array.
00320   * \item[minirn, minicn]
00321   *   Integer: In the event of a
00322   *   successful return (#iflag# >= 0 or #iflag# = -14) give the minimum size
00323   *   of #irn# and #a#/#icn# respectively which would enable a successful run
00324   *   on an identical matrix. On an exit with #iflag# equal to -5, #minicn#
00325   *   gives the minimum value of #icn# for success on subsequent runs on
00326   *   an identical matrix. in the event of failure with #iflag# = -6, -4,
00327   *   -3, -2, or -1, then #minicn# and #minirn# give the minimum value of
00328   *   #licn# and #lirn# respectively which would be required for a
00329   *   successful decomposition up to the point at which the failure
00330   *   occurred.
00331   * \item[irank]
00332   *   Integer:  Gives an upper bound on the rank of the matrix.
00333   * \item[abort1]
00334   *   Logical:  Default value #true#.  If #abort1# is
00335   *   set to #false# then ma28a/ad will decompose structurally singular
00336   *   matrices (including rectangular ones).
00337   * \item[abort2]
00338   *   Logical:   Default value #true#.  If #abort2# is
00339   *   set to #false# then ma28a/ad will decompose numerically singular
00340   *   matrices.
00341   * \item[idisp]
00342   *   Integer[2]:  On output from ma28a/ad, the
00343   *   indices of the diagonal blocks of the factors lie in positions
00344   *   idisp(1) to idisp(2) of #a#/#icn#. This array must be preserved
00345   *   between a call to ma28a/ad and subsequent calls to ma28b/bd,
00346   *   ma28c/cd or ma28i/id.
00347   * \item[tol]
00348   *   Double Precision:  If it is set to a positive value, then any
00349   *   non-zero whose modulus is less than #tol# will be dropped from the
00350   *   factorization.  The factorization will then require less storage
00351   *   but will be inaccurate.  After a run of ma28a/ad with #tol# positive
00352   *   it is not possible to use ma28b/bd and the user is recommended to
00353   *   use ma28i/id to obtain the solution.  The default value for #tol# is
00354   *   0.0.
00355   * \item[themax]
00356   *   Double Precision:  On exit from ma28a/ad, it will hold the
00357   *   largest entry of the original matrix.
00358   * \item[big]
00359   *   Double Precision:  If #lbig# has been set to #true#, #big# will hold
00360   *   the largest entry encountered during the factorization by ma28a/ad
00361   *   or ma28b/bd.
00362   * \item[dxmax]
00363   *   Double Precision:  On exit from ma28i/id, #dxmax# will be set to
00364   *   the largest component of the solution.
00365   * \item[errmax]
00366   *   Double Precision:  On exit from ma28i/id, If #maxit# is
00367   *   positive, #errmax# will be set to the largest component in the
00368   *   estimate of the error.
00369   * \item[dres]
00370   *   Double Precision:  On exit from ma28i/id, if #maxit# is positive,
00371   *   #dres# will be set to the largest component of the residual.
00372   * \item[cgce]
00373   *   Double Precision:  It is used by ma28i/id to check the
00374   *   convergence rate.  if the ratio of successive corrections is
00375   *   not less than #cgce# then we terminate since the convergence
00376   *   rate is adjudged too slow.
00377   * \item[ndrop]
00378   *   Integer:  If #tol# has been set positive, on exit
00379   *   from ma28a/ad, #ndrop# will hold the number of entries dropped from
00380   *   the data structure.
00381   * \item[maxit]
00382   *   Integer:  It is the maximum number of iterations
00383   *   performed by ma28i/id. It has a default value of 16.
00384   * \item[noiter]
00385   *   Integer:  It is set by ma28i/id to the number of
00386   *   iterative refinement iterations actually used.
00387   * \item[nsrch]
00388   *   Integer:  If #nsrch# is set to a value less than #n#,
00389   *   then a different pivot option will be employed by ma28a/ad.  This
00390   *   may result in different fill-in and execution time for ma28a/ad.
00391   *   If #nsrch# is less than or equal to #n#, the workspace array #iw# can be
00392   *   reduced in length.  The default value for nsrch is 32768.
00393   * \item[istart]
00394   *   Integer:  If #istart# is set to a value other than
00395   *   zero, then the user must supply an estimate of the solution to
00396   *   ma28i/id.  The default value for istart is zero.
00397   * \item[lbig]
00398   *   Logical:  If #lbig# is set to #true#, the value of the
00399   *   largest element encountered in the factorization by ma28a/ad or
00400   *   ma28b/bd is returned in #big#.  setting #lbig# to #true#  will
00401   *   increase the time for ma28a/ad marginally and that for ma28b/bd
00402   *   by about 20%.  The default value for #lbig# is #false#.
00403   */
00404 
00405 // @{
00406 //    begin MA28 Common Block Access
00407 
00408 // / Common block with members: #lp#, #mp#, #lblock#, #grow#
00409 static MA28ED_struct &ma28ed_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA28ED,ma28ed);
00410 // / Common block with members: #eps#, #rmin#, #resid#, #irncp#, #icncp#, #minirc#, #minicn#, #irank#, #abort1#, #abort2#
00411 static MA28FD_struct &ma28fd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA28FD,ma28fd);
00412 // / Common block with members: #idisp#
00413 static MA28GD_struct &ma28gd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA28GD,ma28gd);
00414 // / Common block with members: #tol#, #themax#, #big#, #bxmax#, #errmax#, #dres#, #cgce#, #ndrop#, #maxit#, #noiter#, #nsrch#, #istart#, #lbig#
00415 static MA28HD_struct &ma28hd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA28HD,ma28hd);
00416 // / Common block with members: #lp#, #abort1#, #abort2#, #abort3#
00417 static MA30ED_struct &ma30ed_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA30ED,ma30ed);
00418 // / Common block with members: #irncp#, #icncp#, #irank#, #irn#, #icn#
00419 static MA30FD_struct &ma30fd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA30FD,ma30fd);
00420 // / Common block with members: #eps#, #rmin#
00421 static MA30GD_struct &ma30gd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA30GD,ma30gd);
00422 // / Common block with members: #resid#
00423 static MA30HD_struct &ma30hd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA30HD,ma30hd);
00424 // / Common block with members: #tol#, #big#, #ndrop#, #nsrch#, #lbig#
00425 static MA30ID_struct &ma30id_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MA30ID,ma30id);
00426 // / Common block with members: #lp#, #numnz#, #num#, #large#, #abort#
00427 static MC23BD_struct &mc23bd_cb = FORTRAN_COMMMON_BLOCK_NAME_UL(MC23BD,mc23bd);
00428 
00429   // The reason that these are declared static is because I need to
00430   // make sure that these references are initialized before they are
00431   // used in other global defintions.  This means that every translation
00432   // unit will have their own copy of this data.  To reduce this code
00433   // blot you could declair them as pointers then set then using the 
00434   // trick of an initialization class (Myers).
00435 
00436 //    end MA28 Common Block Access
00437 // @}
00438 
00439 // /
00440 /* * Analyze and factor a sparse matrix.
00441   *
00442   * This function analyzes (determines row and column pivots to minimize
00443   * fill-in and result in a better conditioned factorization) then factors
00444   * (calculates upper and lower triangular factors for the determined
00445   * row and column pivots) the permuted system.  The function takes a sparse
00446   * matrix stored in coordinate form ([Aij, i, j] => [#a#, #irn#, #icn#]) and
00447   * factorizes it.  On entry, the first #nz# elements of #a#, #irn#, and
00448   * #icn# hold the matrix elements.  The remaining entires of #a#, #irn#
00449   * , and #icn# hold the fill-in entries after the matrix factorization is
00450   * complete.
00451   *
00452   * The amount of fill-in is influenced by #u#.  A value of #u = 1.0# gives partial
00453   * pivoting where fill-in is sacrificed for the sake of numerical stability and
00454   * #u = 0.0# gives pivoting that strictly minimizes fill-in.
00455   *
00456   * The parameters #ikeep#, #iw#, and #w# are used as workspace but
00457   * #ikeep# contains important information about the factoriation
00458   * on return.
00459   *
00460   * The parameter #iflag# is used return error information about the attempt
00461   * to factorized the system.
00462   *
00463   * @param  n [input] Order of the system being solved.
00464   * @param  nz  [input] Number of non-zeros of the input matrix (#nz >= n#).
00465   *       The ratio #nz/(n*n)# equals the sparsity fraction for the matrix.
00466   * @param  a [input/output] Length = #licn#.  The first #nz# entries hold the
00467   *       non-zero entries of the input matrix on input and
00468   *       the non-zero entries of the factorized matrix on exit.
00469   * @param  licn  [input] length of arrays #a# and #icn#.  This
00470   *         is the total amount of storage advalable for the
00471   *         non-zero entries of the factorization of #a#.
00472   *         therefore #licn# must be greater than #nz#. How
00473   *         much greater depends on the amount of fill-in.
00474   * @param  irn [input/modifed] Length = #ircn#.  The first #nz# entries hold
00475   *       the row indices of the matrix #a# on input.
00476   * @param  ircn  [input] Lenght of irn.
00477   * @param  icn [input/output] Length = #licn#.  Holds column indices of #a# on 
00478   *       entry and the column indices of the reordered
00479   *       #a# on exit.
00480   * @param  u [input] Controls partial pivoting.
00481   *       \begin{description}
00482   *       \item[#* u >= 1.0#]
00483   *         Uses partial pivoting for maximum numerical stability
00484   *         at the expense of some extra fill-in.
00485   *       \item[#* 1.0 < u < 0.0#]
00486   *         Balances numerical stability and fill-in with #u# near 1.0
00487   *         favoring stability and #u# near 0.0 favoring less fill-in.
00488   *       \item[#* u <= 0.0#]
00489   *         Determines row and column pivots to minimize fill-in
00490   *         irrespective of the numerical stability of the 
00491   *         resulting factorization.
00492   *       \end{description}
00493   * @param  ikeep [output] Length = 5 * #n#.  On exist contains information about 
00494   *         the factorization.
00495   *         \begin{description}
00496   *         \item[#* #ikeep(:,1)]
00497   *           Holds the total length of the part of row i in
00498   *           the diagonal block.
00499   *         \item[#* #ikeep(:,2)]
00500   *           Holds the row pivots.  Row #ikeep(i,2)# of the
00501   *           input matrix is the ith row of the pivoted matrix
00502   *           which is factorized.
00503   *         \item[#* #ikeep(:,3)]
00504   *           Holds the column pivots.  Column #ikeep(i,3)# of the
00505   *           input matrix is the ith column of the pivoted matrix
00506   *           which is factorized.
00507   *         \item[#* #ikeep(:,4)]
00508   *           Holds the length of the part of row i in the L
00509   *           part of the L/U decomposition.
00510   *         \item[#* #ikeep(:,5)]
00511   *           Holds the length of the part of row i in the
00512   *           off-diagonal blocks.  If there is only one
00513   *           diagonal block, #ikeep(i,5)# is set to -1.
00514   *         \end{description}
00515   * @param  iw  [] Length = #8*n#.  Integer workspace.
00516   * @param  w [] Length = #n#.  Real workspace.
00517   * @param  iflag [output] Used to return error condtions.
00518   *         \begin{description}
00519   *         \item[#* >= 0#] Success.
00520   *         \item[#* < 0#] Some error has occured.
00521   *         \end{description}
00522   */
00523 inline void ma28ad(const f_int& n, const f_int& nz, f_dbl_prec a[], const f_int& licn
00524   , f_int irn[], const f_int& lirn, f_int icn[], const f_dbl_prec& u, f_int ikeep[], f_int iw[]
00525   , f_dbl_prec w[], f_int* iflag)
00526 { Fortran::FORTRAN_FUNC_CALL_UL(MA28AD,ma28ad) (n,nz,a,licn,irn,lirn,icn,u,ikeep,iw,w,iflag); }
00527 
00528 // /
00529 /* * Factor a sparse matrix using previous analyze pivots.
00530   *
00531   * This function uses the pivots determined from a previous factorization
00532   * to factorize the matrix #a# again.  The assumption is that the 
00533   * sparsity structure of #a# has not changed but only its numerical 
00534   * values.  It is therefore possible that the refactorization may 
00535   * become unstable.
00536   *
00537   * The matrix to be refactorized on order #n# with #nz# non-zero elements
00538   * is input in coordinate format in #a#, #ivect# and #jvect#. 
00539   *
00540   * Information about the factorization is contained in the #icn# and
00541   * #ikeep# arrays returned from \Ref{ma28ad}.
00542   *
00543   * @param  n [input] Order of the system being solved.
00544   * @param  nz  [input] Number of non-zeros of the input matrix
00545   * @param  a [input/output] Length = #licn#.  The first #nz# entries hold the
00546   *       non-zero entries of the input matrix on input and
00547   *       the non-zero entries of the factorized matrix on exit.
00548   * @param  licn  [input] length of arrays #a# and #icn#.  This
00549   *         is the total amount of storage avalable for the
00550   *         non-zero entries of the factorization of #a#.
00551   *         therefore #licn# must be greater than #nz#. How
00552   *         much greater depends on the amount of fill-in.
00553   * @param  icn [input] Length = #licn#.  Same array output from #ma28ad#.
00554   *       It contains information about the analyze step.
00555   * @param  ikeep [input] Length = 5 * #n#.  Same array output form #ma28ad#.
00556   *         It contains information about the analyze step.
00557   * @param  iw  [] Length = #8*n#.  Integer workspace.
00558   * @param  w [] Length = #n#.  Real workspace.
00559   * @param  iflag [output] Used to return error condtions.
00560   *         \begin{description}
00561   *         \item[#* >= 0#] Success
00562   *         \item[#* < 0#] Some error has occured.
00563   *         \end{description}
00564   */
00565 inline void ma28bd(const f_int& n, const f_int& nz, f_dbl_prec a[], const f_int& licn
00566   , const f_int ivect[], const f_int jvect[], const f_int icn[], const f_int ikeep[], f_int iw[]
00567   , f_dbl_prec w[], f_int* iflag)
00568 { Fortran::FORTRAN_FUNC_CALL_UL(MA28BD,ma28bd) (n,nz,a,licn,ivect,jvect,icn,ikeep,iw,w,iflag);  }
00569 
00570 // /
00571 /* * Solve for a rhs using a factorized matrix.
00572   *
00573   * This function solves for a rhs given a matrix factorized by
00574   * #ma28ad# or #ma28bd#.  The right hand side (rhs) is passed
00575   * in in #rhs# and the solution is return in #rhs#.  The 
00576   * factorized matrix is passed in in #a#, #icn#, and #ikeep#
00577   * which were set by #ma28ad# and/or #ma28bd#. The 
00578   *
00579   * The matrix or its transpose can be solved for by selecting
00580   * #mtype == 1# or #mtype != 1# respectively.
00581   *
00582   * @param  n [input] Order of the system being solved.
00583   * @param  a [input] Length = #licn#.  Contains the non-zero
00584   *       elements of the factorized matrix.
00585   * @param  licn  [input] length of arrays #a# and #icn#.  This
00586   *         is the total amount of storage avalable for the
00587   *         non-zero entries of the factorization of #a#.
00588   *         therefore #licn# must be greater than #nz#. How
00589   *         much greater depends on the amount of fill-in.
00590   * @param  icn [input] Length = #licn#.  Same array output from #ma28ad#.
00591   *       It contains information about the analyze step.
00592   * @param  ikeep [input] Length = 5 * #n#.  Same array output form #ma28ad#.
00593   *         It contains information about the analyze step.
00594   * @param  w [] Length = #n#.  Real workspace.
00595   * @param  mtype [input] Instructs to solve using the matrix or its transpoze.
00596   *         \begin{description}
00597   *         \item[#* mtype == 1#] Solve using the non-transposed matrix.
00598   *         \item[#* mtype != 1#] Solve using the transposed matrix.
00599   *         \end{description} 
00600   */
00601 inline void ma28cd(const f_int& n, const f_dbl_prec a[], const f_int& licn, const f_int icn[]
00602   , const f_int ikeep[], f_dbl_prec rhs[], f_dbl_prec w[], const f_int& mtype)
00603 { Fortran::FORTRAN_FUNC_CALL_UL(MA28CD,ma28cd) (n,a,licn,icn,ikeep,rhs,w,mtype);  }
00604 
00605 //    end MA28 C++ Declarations
00606 // @}
00607 
00608 } // end namespace MA28_CDecl
00609 
00610 #endif // MA28_CPPDECL_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends