00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "TQpProbDens.h"
00044 #include "TMatrixD.h"
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 ClassImp(TQpProbDens)
00055
00056
00057 TQpProbDens::TQpProbDens(Int_t nx,Int_t my,Int_t mz) :
00058 TQpProbBase(nx,my,mz)
00059 {
00060
00061
00062
00063 R__ASSERT(nx-my-mz > 0);
00064 }
00065
00066
00067
00068 TQpProbDens::TQpProbDens(const TQpProbDens &another) : TQpProbBase(another)
00069 {
00070
00071
00072 *this = another;
00073 }
00074
00075
00076
00077 TQpDataBase *TQpProbDens::MakeData(Double_t *c,
00078 Double_t *Q,
00079 Double_t *xlo,Bool_t *ixlo,
00080 Double_t *xup,Bool_t *ixup,
00081 Double_t *A, Double_t *bA,
00082 Double_t *C,
00083 Double_t *clo,Bool_t *iclo,
00084 Double_t *cup,Bool_t *icup)
00085 {
00086
00087
00088 TVectorD vc ; vc .Use(fNx,c);
00089 TMatrixDSym mQ ; mQ .Use(fNx,Q);
00090 TVectorD vxlo; vxlo.Use(fNx,xlo);
00091 TVectorD vxup; vxup.Use(fNx,xup);
00092 TMatrixD mA ;
00093 TVectorD vbA ;
00094 if (fMy > 0) {
00095 mA .Use(fMy,fNx,A);
00096 vbA .Use(fMy,bA);
00097 }
00098 TMatrixD mC ;
00099 TVectorD vclo;
00100 TVectorD vcup;
00101 if (fMz > 0) {
00102 mC .Use(fMz,fNx,C);
00103 vclo.Use(fMz,clo);
00104 vcup.Use(fMz,cup);
00105 }
00106
00107 TVectorD vixlo(fNx);
00108 TVectorD vixup(fNx);
00109 for (Int_t ix = 0; ix < fNx; ix++) {
00110 vixlo[ix] = (ixlo[ix]) ? 1.0 : 0.0;
00111 vixup[ix] = (ixup[ix]) ? 1.0 : 0.0;
00112 }
00113
00114 TVectorD viclo(fMz);
00115 TVectorD vicup(fMz);
00116 for (Int_t ic = 0; ic < fMz; ic++) {
00117 viclo[ic] = (iclo[ic]) ? 1.0 : 0.0;
00118 vicup[ic] = (icup[ic]) ? 1.0 : 0.0;
00119 }
00120
00121 TQpDataDens *data = new TQpDataDens(vc,mQ,vxlo,vixlo,vxup,vixup,mA,vbA,mC,vclo,
00122 viclo,vcup,vicup);
00123
00124 return data;
00125 }
00126
00127
00128
00129 TQpDataBase *TQpProbDens::MakeData(TVectorD &c,
00130 TMatrixDBase &Q_in,
00131 TVectorD &xlo, TVectorD &ixlo,
00132 TVectorD &xup, TVectorD &ixup,
00133 TMatrixDBase &A_in,TVectorD &bA,
00134 TMatrixDBase &C_in,
00135 TVectorD &clo, TVectorD &iclo,
00136 TVectorD &cup, TVectorD &icup)
00137 {
00138
00139
00140 TMatrixDSym &mQ = (TMatrixDSym &) Q_in;
00141 TMatrixD &mA = (TMatrixD &) A_in;
00142 TMatrixD &mC = (TMatrixD &) C_in;
00143
00144 R__ASSERT(mQ.GetNrows() == fNx && mQ.GetNcols() == fNx);
00145 if (fMy > 0) R__ASSERT(mA.GetNrows() == fMy && mA.GetNcols() == fNx);
00146 else R__ASSERT(mA.GetNrows() == fMy);
00147 if (fMz > 0) R__ASSERT(mC.GetNrows() == fMz && mC.GetNcols() == fNx);
00148 else R__ASSERT(mC.GetNrows() == fMz);
00149
00150 R__ASSERT(c.GetNrows() == fNx);
00151 R__ASSERT(xlo.GetNrows() == fNx);
00152 R__ASSERT(ixlo.GetNrows() == fNx);
00153 R__ASSERT(xup.GetNrows() == fNx);
00154 R__ASSERT(ixup.GetNrows() == fNx);
00155
00156 R__ASSERT(bA.GetNrows() == fMy);
00157 R__ASSERT(clo.GetNrows() == fMz);
00158 R__ASSERT(iclo.GetNrows() == fMz);
00159 R__ASSERT(cup.GetNrows() == fMz);
00160 R__ASSERT(icup.GetNrows() == fMz);
00161
00162 TQpDataDens *data = new TQpDataDens(c,mQ,xlo,ixlo,xup,ixup,mA,bA,mC,clo,iclo,cup,icup);
00163
00164 return data;
00165 }
00166
00167
00168
00169 TQpResidual* TQpProbDens::MakeResiduals(const TQpDataBase *data_in)
00170 {
00171
00172
00173 TQpDataDens *data = (TQpDataDens *) data_in;
00174 return new TQpResidual(fNx,fMy,fMz,data->fXloIndex,data->fXupIndex,data->fCloIndex,data->fCupIndex);
00175 }
00176
00177
00178
00179 TQpVar* TQpProbDens::MakeVariables(const TQpDataBase *data_in)
00180 {
00181
00182
00183 TQpDataDens *data = (TQpDataDens *) data_in;
00184
00185 return new TQpVar(fNx,fMy,fMz,data->fXloIndex,data->fXupIndex,data->fCloIndex,data->fCupIndex);
00186 }
00187
00188
00189
00190 TQpLinSolverBase* TQpProbDens::MakeLinSys(const TQpDataBase *data_in)
00191 {
00192
00193
00194 TQpDataDens *data = (TQpDataDens *) data_in;
00195 return new TQpLinSolverDens(this,data);
00196 }
00197
00198
00199
00200 void TQpProbDens::JoinRHS(TVectorD &rhs,TVectorD &rhs1_in,TVectorD &rhs2_in,TVectorD &rhs3_in)
00201 {
00202
00203
00204
00205
00206
00207
00208 rhs.SetSub(0,rhs1_in);
00209 if (fMy > 0) rhs.SetSub(fNx, rhs2_in);
00210 if (fMz > 0) rhs.SetSub(fNx+fMy,rhs3_in);
00211 }
00212
00213
00214
00215 void TQpProbDens::SeparateVars(TVectorD &x_in,TVectorD &y_in,TVectorD &z_in,TVectorD &vars_in)
00216 {
00217
00218
00219
00220
00221
00222
00223 x_in = vars_in.GetSub(0,fNx-1);
00224 if (fMy > 0) y_in = vars_in.GetSub(fNx, fNx+fMy-1);
00225 if (fMz > 0) z_in = vars_in.GetSub(fNx+fMy,fNx+fMy+fMz-1);
00226 }
00227
00228
00229
00230 void TQpProbDens::MakeRandomData(TQpDataDens *&data,TQpVar *&soln,Int_t ,Int_t ,Int_t )
00231 {
00232
00233
00234 data = new TQpDataDens(fNx,fMy,fMz);
00235 soln = this->MakeVariables(data);
00236 data->DataRandom(soln->fX,soln->fY,soln->fZ,soln->fS);
00237 }
00238
00239
00240
00241 TQpProbDens &TQpProbDens::operator=(const TQpProbDens &source)
00242 {
00243
00244
00245 if (this != &source) {
00246 TQpProbBase::operator=(source);
00247 }
00248 return *this;
00249 }