Changeset 23b4ae2 in tol


Ignore:
Timestamp:
Feb 3, 2009, 9:14:19 AM (18 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
84e50984
Parents:
1d83266
Message:

New argument forceNaturalOrder in CholeskiFactor
CHOLMOD API for choleski solve treat P and Pt as products instaed of linear sistem solvers

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@559 be4846ae-f051-483a-831f-b64e7d08527f

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol/btol/vmatrix_type/vmat_chol.cpp

    r1d83266 r23b4ae2  
    109109////////////////////////////////////////////////////////////////////////////////
    110110  int BVMat::CholeskiFactor(const BVMat& X_, BVMat& L_,
    111                             ECholFacOri ori,  bool checkSymmetric)
     111                            ECholFacOri ori, 
     112                            bool checkSymmetric,
     113                            bool forceNaturalOrder)
    112114////////////////////////////////////////////////////////////////////////////////
    113115{
    114116  if(!X_.CheckDefined("CholeskiFactor")) { return(-1); }
    115 
     117  bool old_force = force_natural_order(forceNaturalOrder);
    116118  BVMat* X__, *L__;
    117   convertIfNeeded_cRt2cRs(X_,L_,X__,L__,"CholeskiSolve");
     119  convertIfNeeded_cRt2cRs(X_,L_,X__,L__,"CholeskiFactor");
    118120  BVMat &X = *X__,  &L = *L__;
    119121  int r = X.Rows();
    120122  int c = X.Columns();
     123  int result = 0;
    121124  if((ori==ECFO_X)&&(r!=c))
    122125  {
    123126    err_cannot_apply("CholeskiFactor",I2("non square","no cuadrada"),X);
    124     return(-1);
     127    result = -1;
    125128  }
    126129  else if((ori==ECFO_XtX)&&(r<c))
    … …  
    129132      I2("row deficient to build",
    130133         "no tiene bastantes filas para construir")+ " S = X' X",X);
    131     return(-2);
     134    result =-2;
    132135  }
    133136  else if((ori==ECFO_XXt)&&(r>c))
    … …  
    136139      I2("column deficient to build",
    137140         "no tiene bastantes columnas para construir")+" S = X X'",X);
    138     return(-3);
    139   }
    140   if(ori!=ECFO_X) { checkSymmetric=false; }
    141   bool isNotSymm = (checkSymmetric)?!X.IsSymmetric():false;
    142   bool isOk = true;
    143   bool isNotPosDef=false;
    144   if(isNotSymm)
    145   {
    146     err_cannot_apply("CholeskiFactor",
    147       I2("non symmetric","no simetrica"),X);
    148     return(-4);
    149   }
    150   if(checkSymmetric && !isNotSymm) { ((BVMat&)X).CompactSymmetric(false); }
    151   const StrCholFac* cholFac = FindCholFac(X.code_, ori);
    152   if(cholFac)
    153   {
    154     (*cholFac->fun)(X, L, isOk, isNotPosDef);
     141    result =-3;
    155142  }
    156143  else
    157144  {
    158     err_invalid_subtype("CholeskiFactor",X);
    159     isOk = false;
    160   }
    161   if(!isOk)
    162   {
    163     if(isNotPosDef)
     145    if(ori!=ECFO_X) { checkSymmetric=false; }
     146    bool isNotSymm = (checkSymmetric)?!X.IsSymmetric():false;
     147    bool isOk = true;
     148    bool isNotPosDef=false;
     149    if(isNotSymm)
    164150    {
    165151      err_cannot_apply("CholeskiFactor",
    166                        I2("non positive definite",
    167                           "no definida positiva"),X);
    168       L.Delete();
    169       return(-5);
     152        I2("non symmetric","no simetrica"),X);
     153      result = -4;
    170154    }
    171     return(-6);
    172   }   
     155    else
     156    {
     157      if(checkSymmetric && !isNotSymm) { ((BVMat&)X).CompactSymmetric(false); }
     158      const StrCholFac* cholFac = FindCholFac(X.code_, ori);
     159      if(cholFac)
     160      {
     161        (*cholFac->fun)(X, L, isOk, isNotPosDef);
     162      }
     163      else
     164      {
     165        err_invalid_subtype("CholeskiFactor",X);
     166        isOk = false;
     167      }
     168      if(!isOk)
     169      {
     170        if(isNotPosDef)
     171        {
     172          err_cannot_apply("CholeskiFactor",
     173                           I2("non positive definite",
     174                              "no definida positiva"),X);
     175          L.Delete();
     176          result = -5;
     177        }
     178        else
     179        {
     180          result = -6;
     181        }
     182      }
     183    }
     184  }
     185  force_natural_order(old_force);
    173186  return(0);
    174187}
    … …  
    176189////////////////////////////////////////////////////////////////////////////////
    177190  int BVMat::CholeskiFactor(const BVMat& X, BVMat& L,
    178                             const BText& oriName, bool checkSymmetric)
     191                            const BText& oriName,
     192                            bool checkSymmetric,
     193                            bool forceNaturalOrder)
    179194////////////////////////////////////////////////////////////////////////////////
    180195{
    … …  
    190205    return(-3);
    191206  }
    192   return(CholeskiFactor(X,L,ori,checkSymmetric));
    193 }
    194  
    195 ////////////////////////////////////////////////////////////////////////////////
    196   BVMat BVMat::CholeskiFactor(ECholFacOri ori, bool checkSymmetric)
     207  return(CholeskiFactor(X,L,ori,checkSymmetric,forceNaturalOrder));
     208}
     209 
     210////////////////////////////////////////////////////////////////////////////////
     211  BVMat BVMat::CholeskiFactor(ECholFacOri ori,
     212                              bool checkSymmetric,
     213                              bool forceNaturalOrder)
    197214//Matrix algebra operator
    198215////////////////////////////////////////////////////////////////////////////////
    199216{
    200217  BVMat aux;
    201   CholeskiFactor(*this, aux, ori, checkSymmetric);
     218  CholeskiFactor(*this, aux, ori, checkSymmetric, forceNaturalOrder);
    202219  return(aux);
    203220};
    204221
    205222////////////////////////////////////////////////////////////////////////////////
    206   BVMat BVMat::CholeskiFactor(const BText& oriName, bool checkSymmetric)
     223  BVMat BVMat::CholeskiFactor(const BText& oriName,
     224                              bool checkSymmetric,
     225                              bool forceNaturalOrder)
    207226//Matrix algebra operator
    208227////////////////////////////////////////////////////////////////////////////////
    209228{
    210229  BVMat aux;
    211   CholeskiFactor(*this, aux, oriName, checkSymmetric);
     230  CholeskiFactor(*this, aux, oriName, checkSymmetric,forceNaturalOrder);
    212231  return(aux);
    213232}
    … …  
    357376  x.s_.blasRdense_ = cholmod_solve
    358377  (
     378    CHOLMOD_Pt,
     379    L.s_.chlmRfactor_,
     380    b.s_.blasRdense_,
     381    common_
     382  );
     383  return(0);   
     384};
     385
     386////////////////////////////////////////////////////////////////////////////////
     387  int BVMat::cRs_bRd_cholSolPt(const BVMat& L, const BVMat& b, BVMat& x)
     388////////////////////////////////////////////////////////////////////////////////
     389{
     390  x.code_ = ESC_blasRdense;
     391  x.s_.blasRdense_ = cholmod_solve
     392  (
    359393    CHOLMOD_P,
    360394    L.s_.chlmRfactor_,
    … …  
    366400
    367401////////////////////////////////////////////////////////////////////////////////
    368   int BVMat::cRs_bRd_cholSolPt(const BVMat& L, const BVMat& b, BVMat& x)
    369 ////////////////////////////////////////////////////////////////////////////////
    370 {
    371   x.code_ = ESC_blasRdense;
    372   x.s_.blasRdense_ = cholmod_solve
    373   (
    374     CHOLMOD_Pt,
    375     L.s_.chlmRfactor_,
    376     b.s_.blasRdense_,
    377     common_
    378   );
    379   return(0);   
    380 };
    381 
    382 ////////////////////////////////////////////////////////////////////////////////
    383402  int BVMat::cRs_cRs_cholSolPtLLtP(const BVMat& factor, const BVMat& b, BVMat& x)
    384403////////////////////////////////////////////////////////////////////////////////
    … …  
    440459  cholmod_sparse* aux = cholmod_spsolve
    441460  (
    442     CHOLMOD_Pt,
     461    CHOLMOD_P,
    443462    factor.s_.chlmRfactor_,
    444463    b_,
    445464    common_
    446465  );
     466  if(b_ != b.s_.chlmRsparse_)
     467  {
     468    cholmod_free_sparse(&b_, common_);
     469  }
    447470  x.s_.chlmRsparse_= cholmod_spsolve
    448471  (
    … …  
    453476  );
    454477  cholmod_free_sparse(&aux, common_);
    455   if(b_ != b.s_.chlmRsparse_)
    456   {
    457     cholmod_free_sparse(&b_, common_);
    458   }
    459478  return(0);   
    460479};
    … …  
    477496    common_
    478497  );
     498  if(b_ != b.s_.chlmRsparse_)
     499  {
     500    cholmod_free_sparse(&b_, common_);
     501  }
    479502  x.s_.chlmRsparse_= cholmod_spsolve
    480503  (
    481     CHOLMOD_P,
     504    CHOLMOD_Pt,
    482505    factor.s_.chlmRfactor_,
    483506    aux,
    … …  
    485508  );
    486509  cholmod_free_sparse(&aux, common_);
    487   if(b_ != b.s_.chlmRsparse_)
    488   {
    489     cholmod_free_sparse(&b_, common_);
    490   }
    491510  return(0);   
    492511};
    … …  
    552571  x.s_.chlmRsparse_ = cholmod_spsolve
    553572  (
    554     CHOLMOD_P,
     573    CHOLMOD_Pt,
    555574    L.s_.chlmRfactor_,
    556575    b_,
    … …  
    576595  x.s_.chlmRsparse_ = cholmod_spsolve
    577596  (
    578     CHOLMOD_Pt,
     597    CHOLMOD_P,
    579598    L.s_.chlmRfactor_,
    580599    b_,
Note: See TracChangeset for help on using the changeset viewer.