00001 #ifndef PRIMALDUALSYSTEM
00002 #define PRIMALDUALSYSTEM
00003
00004 #include "IotrRefCount.hh"
00005 #include "IotrHandle.hh"
00006 #include "LinearOperator.hh"
00007 #include "IotrUtilities.hh"
00008
00009 class GenMatrix;
00010 class SymMatrix;
00011 typedef Handle<SymMatrix> SymMatrixHandle;
00012
00013 class IotrVector;
00014 typedef Handle<IotrVector> IotrVectorHandle;
00015
00016 class GenMatrix;
00017 typedef Handle<GenMatrix> GenMatrixHandle;
00018
00019 class LinearSymSolver;
00020 typedef Handle<LinearSymSolver> LinearSymSolverHandle;
00021
00022 #include "ElementGather.hh"
00023
00029 class PrimalDualSystem : public IotrRefCount {
00030 friend class PrimalDualSystemTester;
00031 protected:
00032 IotrVectorHandle mDx;
00033 IotrVectorHandle mOmega;
00034 IotrVectorHandle mDgamma;
00035 IotrVectorHandle mDphi;
00036 IotrVectorHandle mDlambda;
00037 IotrVectorHandle mDpi;
00038 IotrVectorHandle mRc;
00039
00040 ElementGatherHandle mIclow;
00041 ElementGatherHandle mIcupp;
00042 ElementGatherHandle mIeq;
00043 ElementGatherHandle mIxlow;
00044 ElementGatherHandle mIxupp;
00045
00046 LinearOperatorHandle mJmult;
00047
00048 int mState;
00049 public:
00050 PrimalDualSystem( ElementGather * iclow, ElementGather * icupp,
00051 ElementGather * ieq,
00052 ElementGather * ixlow, ElementGather * ixupp );
00053
00054
00055 virtual void setDiagonals( IotrVector & Dx, IotrVector & Dy,
00056 IotrVector & Dgamma, IotrVector & Dphi,
00057 IotrVector & Dlambda, IotrVector & Dpi );
00058
00059 virtual void assemble( GenMatrix & J, SymMatrix & H,
00060 IotrVector & Dx, IotrVector & omega ) = 0;
00061
00062 virtual void factor( IotrVector & Dx, IotrVector & Dy,
00063 IotrVector & Dgamma, IotrVector & Dphi,
00064 IotrVector & Dlambda, IotrVector & Dpi,
00065 LinearOperator * Jmult,
00066 GenMatrix & J, SymMatrix & H );
00067
00068 virtual void matrixChanged() = 0;
00069
00070 virtual void newtonStep( IotrVector & xsol,
00071 IotrVector & gammaSol, IotrVector & phiSol,
00072 IotrVector & ysol,
00073 IotrVector & lambdaSol, IotrVector & piSol );
00074
00075 virtual void basicSolve( IotrVector & xsol,
00076 IotrVector & rc ) = 0;
00077
00078 virtual int isConvex() = 0;
00079 virtual void inspect();
00080 virtual int isValid();
00081
00082 int isSingular()
00083 {
00084 IotrAssert( mState == factored || mState == singular );
00085
00086 return mState == singular;
00087 }
00088
00089 ElementGather & ixlow() { return *mIxlow; }
00090 ElementGather & ixupp() { return *mIxupp; }
00091 ElementGather & ieq() { return *mIeq; }
00092 ElementGather & iclow() { return *mIclow; }
00093 ElementGather & icupp() { return *mIcupp; }
00094
00095 int mc() { return mIeq->nx(); }
00096 int nx() { return mIxlow->nx(); }
00097 int my() { return mIeq->my(); }
00098 int mclow() { return mIclow->my(); }
00099 int mcupp() { return mIcupp->my(); }
00100 int nxlow() { return mIxlow->my(); }
00101 int nxupp() { return mIxupp->my(); }
00102
00103 IotrVector & omega() { return *mOmega; }
00104 IotrVector & Dx() { return *mDx; }
00105
00106 enum {initialState = 0, diagonalsSet, assembled, factored, singular};
00107
00108 int state() { return mState; }
00109 };
00110
00111 typedef Handle<PrimalDualSystem> PrimalDualSystemHandle;
00112
00113 #endif
