1 #include "LinearEquation.h"
2
3 using NEWMAT::ColumnVector;
4 using NEWMAT::Matrix;
5 using NEWMAT::SymmetricMatrix;
6
7 //------------------------------------------------------------------------
8 // P.J. Williams
9 // Sandia National Laboratories
10 // pwillia@sandia.gov
11 // Last modified 03/20/2003
12 //------------------------------------------------------------------------
13
14 namespace OPTPP {
15
16 // Constructors
LinearEquation()17 LinearEquation::LinearEquation():
18 LinearConstraint(), b_(0), ctype_(0){;}
19
LinearEquation(const Matrix & A,const ColumnVector & rhs)20 LinearEquation::LinearEquation(const Matrix& A, const ColumnVector& rhs):
21 LinearConstraint(A, rhs), b_(rhs), ctype_(A.Nrows())
22 {ctype_.ReSize(numOfCons_); ctype_ = Leqn;}
23
24 // Functions for computing various quantities
evalAx(const ColumnVector & xc)25 ColumnVector LinearEquation::evalAx(const ColumnVector& xc) const
26 {
27 int i, index;
28 ColumnVector Ax(numOfCons_);
29 Matrix tmp(numOfCons_, numOfVars_);
30 for( i = 1; i <= numOfCons_; i++){
31 index = constraintMappingIndices_[i-1];
32 tmp.Row(i) = A_.Row(index);
33 }
34 Ax = tmp*xc;
35 return Ax;
36 }
37
evalCFGH(const ColumnVector & xc)38 void LinearEquation::evalCFGH(const ColumnVector & xc) const
39 {
40 return;
41 }
42
evalResidual(const ColumnVector & xc)43 ColumnVector LinearEquation::evalResidual(const ColumnVector& xc) const
44 {
45 int i, index;
46 cvalue_ = A_*xc;
47 ColumnVector Ax = evalAx(xc);
48 ColumnVector residual(numOfCons_);
49
50 for( i = 1; i <= numOfCons_; i++){
51 index = constraintMappingIndices_[i-1];
52 residual(i) = Ax(i) - b_(index);
53 }
54 return residual;
55 }
56
evalGradient(const ColumnVector & xc)57 Matrix LinearEquation::evalGradient(const ColumnVector& xc) const
58 {
59 int i, index;
60 Matrix tmp(numOfCons_, numOfVars_);
61
62 for( i = 1; i <= numOfCons_; i++){
63 index = constraintMappingIndices_[i-1];
64 tmp.Row(i) = A_.Row(index);
65 }
66 return tmp.t();
67 }
68
69
amIFeasible(const ColumnVector & xc,double epsilon)70 bool LinearEquation::amIFeasible(const ColumnVector & xc, double epsilon) const
71 {
72 int i;
73 bool feasible = true;
74 ColumnVector residual = evalResidual(xc);
75 for(i = 1; i <= numOfCons_; i++){
76 if( (residual(i) > epsilon) || (residual(i) < -epsilon) ){
77 feasible = false;
78 break;
79 }
80 }
81 return feasible;
82 }
83
84 } // namespace OPTPP
85