- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/modules/solvLin.cc File Reference

`#include <iostream>`
`#include <fstream>`
`#include <math.h>`
`#include <stdlib.h>`
`#include "../def.h"`
`#include "../model.h"`
`#include "../nr.h"`

Go to the source code of this file.

## Defines
#define | PRINTBIGSYS
## Functions
int | lsei_ (double *w, long &mdw, long &me, long &ma, long &mg, long &n, double *prgopt, double *x, double *rnorme, double *rnorml, long *mode, double *ws, long *ip)
int | e04ncf_ (long &M, long &N, long &NCLIN, long &LDC, long &LDA, double *C, double *BL, double *BU, double *CVEC, long *ISTATE, long *KX, double *X, double *A, double *B, long &ITER, double &OBJ, double *CLAMBDA, long *IWORK, long &LIWORK, double *WORK, long &LWORK, long &IFAIL)
int | setopte04ncf_ ()
void | condense (Glob *globs, GlobExp *ex)
| Condenses each experiment.
void | decondense (Glob *globs, GlobExp *ex, double *dX)
| Decondenses each experiment.
void | solvLin (Glob *globs, GlobExp *ex, int computeCovar)
void | myerr_ (int *messg, int &nmessg)
void | xerr_ (char *msg, int &nerr, int &iert, int &ni, int &i1, int &i2, int &nr, double &r1, double &r2)

## Define Documentation

#define PRINTBIGSYS

| Definition at line 14 of file solvLin.cc.

## Function Documentation

void condense | ( | Glob * | globs,
|  | GlobExp * | ex
| ) |

| Condenses each experiment. Definition at line 40 of file solvLin.cc. References dmatrix(), Glob::elastic, GlobExp::me, GlobExp::mg, GlobExp::nobs, Glob::npar, GlobExp::nPoints, Glob::nrExp, and GlobExp::nvar. Referenced by solvLin().
    00042 {
    00043   long nExp,msP,i,j,k,l;
    00044   long nMeas,nvar,nP;
    00045   long me,mg,nobs,nms;
    00046
    00047   //weakens continuty constraints
    00048   double elastic=globs->elastic;
    00049   double **Ea,**Ee,**Eg; //save the previous E-matrices
    00050
    00051   nP=globs->npar;
    00052   for(nExp=1;nExp<=globs->nrExp;++nExp)
    00053     {
    00054       nMeas=ex[nExp].nMeasure;
    00055       nvar=ex[nExp].nvar;
    00056       me=ex[nExp].me;
    00057       mg=ex[nExp].mg;
    00058       nobs=ex[nExp].nobs;
    00059       nms=ex[nExp].nPoints;
    00060       nms--;
    00061       Ea=dmatrix(1,nMeas*nobs,1,nvar);
    00062       Ee=dmatrix(1,me,1,nvar);
    00063       Eg=dmatrix(1,mg,1,nvar);
    00064
    00065       //initialisation of backward recursion
    00066       //------------------------------------
    00067       //least squares:
    00068       for(i=1;i<=nMeas;i++)
    00069         {
    00070           for(j=1;j<=nobs;j++)
    00071             {
    00072               //ua
    00073               ex[nExp].ua[(j-1)*nMeas+i]=ex[nExp].residues[i][j]/ex[nExp].sigma[i][j];
    00074               //Ea
    00075               for(k=1;k<=nvar;k++)
    00076                 {
    00077                   //find out which datapoints are falling in the last MS-interval
    00078                   //the other derivatives are zero
    00079                   if(ex[nExp].xMeasure[i] >= ex[nExp].mesh[nms] &&
    00080                      ex[nExp].xMeasure[i] <= ex[nExp].mesh[nms+1])
    00081                     {
    00082                       ex[nExp].Ea[(j-1)*nMeas+i][k]=ex[nExp].dmds[i][j][k]/ex[nExp].sigma[i][j];
    00083                     }
    00084                   else
    00085                     {
    00086                       ex[nExp].Ea[(j-1)*nMeas+i][k]=0.0;
    00087                     }
    00088                   Ea[(j-1)*nMeas+i][k]=ex[nExp].Ea[(j-1)*nMeas+i][k];
    00089                 }
    00090               //Pa
    00091               for(k=1;k<=nP;k++)
    00092                 {
    00093                   ex[nExp].Pa[(j-1)*nMeas+i][k]=ex[nExp].dmdp[i][j][k]/ex[nExp].sigma[i][j];
    00094                 }
    00095             }
    00096         }// end of least squares init.
    00097
    00098       //equality constraints
    00099       for(i=1;i<=me;i++)
    00100         {
    00101           //ue
    00102           ex[nExp].ue[i]=ex[nExp].r2[i];
    00103           //Ee
    00104           for(j=1;j<=nvar;j++)
    00105             {
    00106               ex[nExp].Ee[i][j]=ex[nExp].dR2ds[i][nms][j];
    00107               Ee[i][j]=ex[nExp].Ee[i][j];
    00108             }
    00109           //Pe
    00110           for(j=1;j<=nP;j++)
    00111             ex[nExp].Pe[i][j]=ex[nExp].dR2dp[i][j];
    00112         }//end of equality constraints
    00113
    00114       //inequality constraints
    00115       for(i=1;i<=mg;i++)
    00116         {
    00117           //ug
    00118           ex[nExp].ug[i]=ex[nExp].r3[i];
    00119           //Ee
    00120           for(j=1;j<=nvar;j++)
    00121             {
    00122               ex[nExp].Eg[i][j]=ex[nExp].dR3ds[i][nms][j];
    00123               Eg[i][j]=ex[nExp].Eg[i][j];
    00124             }
    00125           //Pe
    00126           for(j=1;j<=nP;j++)
    00127             ex[nExp].Pg[i][j]=ex[nExp].dR3dp[i][j];
    00128         }//end of inequality constraints
    00129       //initialisation of backward recursion complete
    00130
    00131
    00132       //backward recursion
    00133       //-----------------
    00134
    00135       //loop over all multiple shooting intervals
    00136       for(msP=nms;msP>=2;msP--)
    00137         {
    00138           //least squares
    00139
    00140           for(i=1;i<=nMeas;i++)
    00141             {
    00142               for(j=1;j<=nobs;j++)
    00143                 {
    00144                   //ua
    00145                   for(k=1;k<=nvar;k++)
    00146                     ex[nExp].ua[(j-1)*nMeas+i]+=ex[nExp].Ea[(j-1)*nMeas+i][k]*ex[nExp].h[msP][k]*elastic;
    00147                   //Ea
    00148                   for(k=1;k<=nvar;k++)
    00149                     {
    00150                       //find out which datapoints are falling in the MS-interval
    00151                       //the other derivatives are zero
    00152                       if(ex[nExp].xMeasure[i] >= ex[nExp].mesh[msP-1] &&
    00153                          ex[nExp].xMeasure[i] <= ex[nExp].mesh[msP])
    00154                         {
    00155                           ex[nExp].Ea[(j-1)*nMeas+i][k]=ex[nExp].dmds[i][j][k]/ex[nExp].sigma[i][j];
    00156                         }
    00157                       else
    00158                         {
    00159                           ex[nExp].Ea[(j-1)*nMeas+i][k]=0.0;
    00160                         }
    00161                       for(l=1;l<=nvar;l++)
    00162                         ex[nExp].Ea[(j-1)*nMeas+i][k]+=Ea[(j-1)*nMeas+i][l]*ex[nExp].dyds[msP][l][k];
    00163                     }
    00164                   //Pa
    00165                   for(k=1;k<=nP;k++)
    00166                     {
    00167                       for(l=1;l<=nvar;l++)
    00168                         {
    00169                           ex[nExp].Pa[(j-1)*nMeas+i][k]+=Ea[(j-1)*nMeas+i][l]*ex[nExp].dydp[msP][l][k];
    00170                         }
    00171                     }
    00172                 }
    00173             } //end of least squares
    00174
    00175           //equality constraints
    00176           for(i=1;i<=me;i++)
    00177             {
    00178               //ue
    00179               for(k=1;k<=nvar;k++)
    00180                 ex[nExp].ue[i]+=Ee[i][k]*ex[nExp].h[msP][k]*elastic;
    00181               //Ee
    00182               for(j=1;j<=nvar;j++)
    00183                 {
    00184                   ex[nExp].Ee[i][j]=ex[nExp].dR2ds[i][msP-1][j];
    00185                     for(l=1;l<=nvar;l++)
    00186                       ex[nExp].Ee[i][j]+=Ee[i][l]*ex[nExp].dyds[msP][l][j];
    00187                 }
    00188               //Pe
    00189               for(j=1;j<=nP;j++)
    00190                 {
    00191                   for(l=1;l<=nvar;l++)
    00192                     ex[nExp].Pe[i][j]+=Ee[i][l]*ex[nExp].dydp[msP][l][j];
    00193                 }
    00194             }//end of equality constraints
    00195
    00196           //inequality constraints
    00197           for(i=1;i<=mg;i++)
    00198             {
    00199               //ue
    00200               for(k=1;k<=nvar;k++)
    00201                 ex[nExp].ug[i]+=Eg[i][k]*ex[nExp].h[msP][k]*elastic;
    00202               //Ee
    00203               for(j=1;j<=nvar;j++)
    00204                 {
    00205                   ex[nExp].Eg[i][j]=ex[nExp].dR3ds[i][msP-1][j];
    00206                     for(l=1;l<=nvar;l++)
    00207                       ex[nExp].Eg[i][j]+=Eg[i][l]*ex[nExp].dyds[msP][l][j];
    00208                 }
    00209               //Pe
    00210               for(j=1;j<=nP;j++)
    00211                 {
    00212                   for(l=1;l<=nvar;l++)
    00213                     ex[nExp].Pg[i][j]+=Eg[i][l]*ex[nExp].dydp[msP][l][j];
    00214                 }
    00215             }//end of inequality constraints
    00216
    00217           //updating new Ea,Ee,Ea
    00218           for(i=1;i<=nvar;i++)
    00219             {
    00220               for(j=1;j<=nMeas*nobs;j++)
    00221                 Ea[j][i]=ex[nExp].Ea[j][i];
    00222               for(j=1;j<=me;j++)
    00223                 Ee[j][i]=ex[nExp].Ee[j][i];
    00224               for(j=1;j<=mg;j++)
    00225                 Eg[j][i]=ex[nExp].Eg[j][i];
    00226             }
    00227
    00228         }//end loop over multiple shooting intervals
    00229
    00230       //clearing memory
    00231       free_dmatrix(Ea,1,nMeas*nobs,1,nvar);
    00232       free_dmatrix(Ee,1,me,1,nvar);
    00233       free_dmatrix(Eg,1,mg,1,nvar);
    00234     }//end loop over all experiments
    00235 }

void decondense | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | double * | dX
| ) |

| Decondenses each experiment. Definition at line 238 of file solvLin.cc. References GlobExp::dS, Glob::elastic, FALSE, Glob::npar, GlobExp::nPoints, Glob::nrExp, GlobExp::nvar, and Glob::y0fix.
    00240 {
    00241   long nExp,msP,i,j,k,l;
    00242   long nvar,npar;
    00243   long nms,ind=0;
    00244
    00245   //weakens continuty constraints
    00246   double elastic=globs->elastic;
    00247
    00248   npar=globs->npar;
    00249   nvar=ex[1].nvar;
    00250
    00251   //initialise
    00252   for(nExp=1;nExp<=globs->nrExp;++nExp)
    00253     {
    00254       for(i=1;i<=ex[nExp].nPoints;i++)
    00255         for(j=1;j<=nvar;j++)
    00256           ex[nExp].dS[i][j]=0.;
    00257       for(i=1;i<=npar;i++)
    00258         ex[nExp].dP[i]=0.;
    00259
    00260       for(i=1;i<=nvar;i++)
    00261         {
    00262           if(globs->y0fix[i]!=FALSE)
    00263             {
    00264               ex[nExp].dS[1][i]=dX[ind];
    00265               ind++;
    00266             }
    00267         }
    00268       for(i=1;i<=npar;i++)
    00269         {
    00270           if(globs->doP[i]=='L')
    00271             {
    00272               ex[nExp].dP[i]=dX[ind];
    00273               ind++;
    00274             }
    00275         }
    00276     }
    00277
    00278   for(i=1;i<=npar;i++)
    00279     {
    00280       if(globs->doP[i]==TRUE)
    00281         {
    00282           for(nExp=1;nExp<=globs->nrExp;++nExp)
    00283             {
    00284               ex[nExp].dP[i]=dX[ind];
    00285             }
    00286           ind++;
    00287         }
    00288     }
    00289
    00290   //end initialise
    00291
    00292   //forward iteration
    00293   for(nExp=1;nExp<=globs->nrExp;++nExp)
    00294     {
    00295       nms=ex[nExp].nPoints;
    00296       nms--;
    00297
    00298       for(msP=2;msP<=nms;msP++)
    00299         {
    00300           for(i=1;i<=nvar;i++)
    00301             {
    00302               ex[nExp].dS[msP][i]=elastic*ex[nExp].h[msP][i];
    00303               for(j=1;j<=nvar;j++)
    00304                 ex[nExp].dS[msP][i]+=ex[nExp].dyds[msP][i][j]*ex[nExp].dS[msP-1][j];
    00305               for(j=1;j<=npar;j++)
    00306                 ex[nExp].dS[msP][i]+=ex[nExp].dydp[msP][i][j]*ex[nExp].dP[j];
    00307             }
    00308         }
    00309
    00310     }//end loop over experiments
    00311 }

int e04ncf_ | ( | long & | M,
|  | long & | N,
|  | long & | NCLIN,
|  | long & | LDC,
|  | long & | LDA,
|  | double * | C,
|  | double * | BL,
|  | double * | BU,
|  | double * | CVEC,
|  | long * | ISTATE,
|  | long * | KX,
|  | double * | X,
|  | double * | A,
|  | double * | B,
|  | long & | ITER,
|  | double & | OBJ,
|  | double * | CLAMBDA,
|  | long * | IWORK,
|  | long & | LIWORK,
|  | double * | WORK,
|  | long & | LWORK,
|  | long & | IFAIL
| ) |

int lsei_ | ( | double * | w,
|  | long & | mdw,
|  | long & | me,
|  | long & | ma,
|  | long & | mg,
|  | long & | n,
|  | double * | prgopt,
|  | double * | x,
|  | double * | rnorme,
|  | double * | rnorml,
|  | long * | mode,
|  | double * | ws,
|  | long * | ip
| ) |

void myerr_ | ( | int * | messg,
|  | int & | nmessg
| ) |

| Definition at line 1034 of file solvLin.cc.
    01035 {
    01036   char *s = (char *) messg;
    01037   cerr << endl;
    01038   for (int i = 0; i < nmessg; i++)
    01039     cerr << s[i];
    01040   cerr << endl;
    01041 }

int setopte04ncf_ | ( |  | ) |

void solvLin | ( | Glob * | globs,
|  | GlobExp * | ex,
|  | int | computeCovar
| ) |

| Definition at line 313 of file solvLin.cc. References condense(), FALSE, Glob::npar, Glob::nrExp, GlobExp::nvar, and Glob::y0fix.
    00314 {
    00315
    00316   long i,j,k,l,xind=0,yind=0,nL;
    00317   long nvarFit=0;  // # der zu fittenden Anfangswerte
    00318   long nLpara=0;   // # lokale Parameter
    00319   long nparL=0;    // nLpara*nrExp
    00320   long npar=globs->npar;
    00321   long nvar=ex[1].nvar;
    00322   long nrExp=globs->nrExp;
    00323   long nparG=0;
    00324   long fitDim=0;
    00325   long aDim=0,eDim=0,gDim=0;
    00326
    00327   //big system (condensed)
    00328   double **Ma,*Ra;
    00329   double **Me,*Re;
    00330   double **Mg,*Rg;
    00331
    00332   //solution
    00333   double *dX;
    00334   //condense experiments
    00335   condense(globs,ex);
    00336
    00337   //determine dimension of big system
    00338   for(i=1;i<=nvar;i++)
    00339     {
    00340       if(globs->y0fix[i]!=FALSE)
    00341         nvarFit++;
    00342     }
    00343   nvarFit*=nrExp;
    00344   for(i=1;i<=npar;i++)
    00345     {
    00346       if(globs->doP[i]=='L')
    00347         nparL++;
    00348       else if(globs->doP[i]!=FALSE)
    00349         nparG++;
    00350     }
    00351   nLpara=nparL;
    00352   nparL*=nrExp;
    00353   fitDim=nvarFit+nparL+nparG;
    00354   for(i=1;i<=nrExp;i++)
    00355     {
    00356       aDim+=ex[i].nMeasure*ex[i].nobs;
    00357       eDim+=ex[i].me;
    00358       gDim+=ex[i].mg;
    00359     }
    00360   if(fitDim==0)
    00361     {
    00362       cerr << "Nothing to fit.\n";
    00363       throw 1;
    00364     }
    00365 #ifdef PRINTDIMENSION
    00366   //print dimension of bigsys;
    00367   *dbg << "Dimension of Big-System\n";
    00368   *dbg << "-----------------------\n";
    00369   *dbg << "Least-squares:\n";
    00370   *dbg << "#Rows : " << aDim  << "  #Columns : " << fitDim << endl;
    00371   *dbg << "Equality constraints:\n";
    00372   *dbg << "#Rows : " << eDim  << "  #Columns : " << fitDim << endl;
    00373   *dbg << "Inequality constraints:\n";
    00374   *dbg << "#Rows : " << gDim  << "  #Columns : " << fitDim;
    00375   *dbg << "\n\n";
    00376   dbg->flush();
    00377 #endif
    00378
    00379
    00380   long aDimold=aDim;
    00381   long eDimold=eDim;
    00382   long fitDimold=fitDim;
    00383
    00384   // Beginn Local_Parameter_Constraints
    00385   if(globs->faktorLexist)
    00386   {
    00387         // Erweiterung der Dimensionen um nLpara*nExp
    00388         aDim=aDim+nparL;
    00389         eDim=eDim+nparL;
    00390         fitDim=fitDim+nparL;
    00391   }
    00392   // Ende Local_Parameter_Constraints
    00393
    00394   //allocate memory
    00395   Ma=dmatrix(0,aDim,1,fitDim);
    00396   Me=dmatrix(0,eDim,1,fitDim);
    00397   Mg=dmatrix(0,gDim,1,fitDim);
    00398   Ra=dvector(0,aDim);
    00399   Re=dvector(0,eDim);
    00400   Rg=dvector(0,gDim);
    00401
    00402   //initialise
    00403   for(i=1;i<=aDim;i++)
    00404     {
    00405       Ra[i]=0.;
    00406       for(j=1;j<=fitDim;j++)
    00407         Ma[i][j]=0.;
    00408     }
    00409    for(i=1;i<=eDim;i++)
    00410     {
    00411       Re[i]=0.;
    00412       for(j=1;j<=fitDim;j++)
    00413         Me[i][j]=0.;
    00414     }
    00415    for(i=1;i<=gDim;i++)
    00416      {
    00417       Rg[i]=0.;
    00418       for(j=1;j<=fitDim;j++)
    00419         Mg[i][j]=0.;
    00420     }
    00421
    00422   //least-squares
    00423   for(i=1;i<=nrExp;i++)
    00424     {
    00425       for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
    00426         Ra[yind+k]=-ex[i].ua[k];
    00427       //initial values
    00428       for(j=1;j<=nvar;j++)
    00429         {
    00430           if(globs->y0fix[j]!=FALSE)
    00431             {
    00432               xind++;
    00433               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
    00434                 Ma[yind+k][xind]=ex[i].Ea[k][j];
    00435             }
    00436
    00437         }
    00438       //parameters
    00439       nL=0;
    00440       for(j=1;j<=npar;j++)
    00441         {
    00442           if(globs->doP[j]=='L')
    00443             {
    00444               xind++;
    00445               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
    00446                 Ma[yind+k][xind]=ex[i].Pa[k][j];
    00447             }
    00448           else if(globs->doP[j]!=FALSE)
    00449             {
    00450               nL++;
    00451               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
    00452                 Ma[yind+k][nvarFit+nparL+nL]=ex[i].Pa[k][j];
    00453             }
    00454         }
    00455       yind=yind+(ex[i].nMeasure*ex[i].nobs);
    00456     }//end loop over experiments
    00457   //-------------> end lest squares
    00458
    00459   //equality constraints
    00460   yind=0;
    00461   xind=0;
    00462   for(i=1;i<=nrExp;i++)
    00463     {
    00464       for(k=1;k<=ex[i].me;k++)
    00465         Re[yind+k]=-ex[i].ue[k];
    00466       //initial values
    00467       for(j=1;j<=nvar;j++)
    00468         {
    00469           if(globs->y0fix[j]!=FALSE)
    00470             {
    00471               xind++;
    00472               for(k=1;k<=ex[i].me;k++)
    00473                 Me[yind+k][xind]=ex[i].Ee[k][xind];
    00474             }
    00475         }
    00476       //parameters
    00477       nL=0;
    00478       for(j=1;j<=npar;j++)
    00479         {
    00480           if(globs->doP[j]=='L')
    00481             {
    00482               xind++;
    00483               for(k=1;k<=ex[i].me;k++)
    00484                 Me[yind+k][xind]=ex[i].Pe[k][j];
    00485             }
    00486           else if(globs->doP[j]!=FALSE)
    00487             {
    00488               nL++;
    00489               for(k=1;k<=ex[i].me;k++)
    00490                 Me[yind+k][nvarFit+nparL+nL]=ex[i].Pe[k][j];
    00491             }
    00492         }
    00493       yind+=ex[i].me;
    00494     }//end loop over experiments
    00495   //-------------> end equality  constraints
    00496
    00497   //inequality constraints
    00498   yind=0;
    00499   xind=0;
    00500   for(i=1;i<=nrExp;i++)
    00501     {
    00502       for(k=1;k<=ex[i].mg;k++)
    00503         Rg[yind+k]=-ex[i].ug[k];
    00504       //initial values
    00505       for(j=1;j<=nvar;j++)
    00506         {
    00507           if(globs->y0fix[j]!=FALSE)
    00508             {
    00509               xind++;
    00510               for(k=1;k<=ex[i].mg;k++)
    00511                 Mg[yind+k][xind]=ex[i].Eg[k][xind];
    00512             }
    00513         }
    00514       //parameters
    00515       nL=0;
    00516       for(j=1;j<=npar;j++)
    00517         {
    00518           if(globs->doP[j]=='L')
    00519             {
    00520               xind++;
    00521               for(k=1;k<=ex[i].mg;k++)
    00522                 Mg[yind+k][xind]=ex[i].Pg[k][j];
    00523             }
    00524           else if(globs->doP[j]!=FALSE)
    00525             {
    00526               nL++;
    00527               for(k=1;k<=ex[i].mg;k++)
    00528                 Mg[yind+k][nvarFit+nparL+nL]=ex[i].Pg[k][j];
    00529             }
    00530         }
    00531       yind+=ex[i].mg;
    00532     }//end loop over experiments
    00533   //-------------> end equality  constraints
    00534
    00535   //Big System Ready to Solve
    00536   //-------------------------
    00537   //
    00538   //  Ma*dX  ~  Ra
    00539   //  Me*dX  =  Re
    00540   //  Mg*dX >=  Rg
    00541   //
    00542   //-------------------------
    00543
    00544   //****************************************************************************
    00545   // Beginn Local_Parameter_Constraints
    00546   // Erweiterung der Systems auf Local_Parameter_Constraints:
    00547   //
    00548   // pLmean und pLsig bestimmen
    00549   if(globs->faktorLexist)
    00550   {
    00551         double *pLmean, *pLsig;
    00552         int index;
    00553         int index_local_parameter;
    00554         pLmean=dvector(1,nLpara);
    00555         pLsig=dvector(1,nLpara);
    00556
    00557 //      for(i=1;i<=nLpara;i++)
    00558 //      {
    00559 //         cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl;
    00560 //         if(globs->faktorL[i]<0)
    00561 //         {
    00562 //            //i;
    00563 //            break;
    00564 //         }
    00565 //      }
    00566 //      cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl;
    00567 //      if(i<=nLpara)
    00568 //      {
    00569 //         cerr<<"number of local parameter constraints too small"<<endl;
    00570 //         throw(1);
    00571 //      }
    00572         for(i=1;i<=nLpara;i++)
    00573         {
    00574                 pLmean[i]=0;
    00575                 index=0;
    00576                 index_local_parameter=0;
    00577                 while(index_local_parameter!=i)
    00578                 {
    00579                         index++;
    00580 //                      cerr<<"doP[index]="<<globs->doP[index]<<" index="<<index<<" i="<<i<<endl;
    00581                         if(globs->doP[index]=='L')
    00582                         {
    00583                                 index_local_parameter++;
    00584 //                              cerr<<"index="<<index<<" i="<<i<<endl;
    00585                         }
    00586
    00587                 }
    00588                 for(j=1;j<=nrExp;j++)
    00589                 {
    00590                         pLmean[i]=pLmean[i]+ex[j].par[index];
    00591 //                      cerr<<"ex[j].par[index]="<<ex[j].par[index]<<endl;
    00592                 }
    00593                 pLmean[i]=pLmean[i]/double(nrExp);
    00594           //    cerr<<"pLmean["<<i<<"]="<<pLmean[i]<<" globs->faktorL["<<i<<"]"<<globs->faktorL[i]<<endl;
    00595         }
    00596         // Berechnung pLsig
    00597         for(i=1;i<=nLpara;i++)
    00598         {
    00599                 pLsig[i]=(globs->faktorL[i]*pLmean[i]-pLmean[i])/(1+globs->faktorL[i]);
    00600                 if(pLsig[i]<0)
    00601                 {
    00602                 cerr<<"ERROR - Negative varianz bei den Lokalen Parameter Constraints..."<<endl;
    00603                 }
    00604         }
    00605         // Ra,Re
    00606         k=1;
    00607         for(i=1;i<=nrExp;i++)
    00608         {
    00609                 for(j=1;j<=nLpara;j++)
    00610                 {
    00611                         index=0;
    00612                         index_local_parameter=0;
    00613                         while(index_local_parameter!=j)
    00614                         {
    00615                                 index++;
    00616                                 if(globs->doP[index]=='L')
    00617                                 {
    00618                                         index_local_parameter++;
    00619                                 }
    00620
    00621                         }
    00622                         Ra[aDimold+k]=(ex[i].par[index]-pLmean[j])/pLsig[j];
    00623                         Re[eDimold+k]=0;
    00624                         double pLsig_real=0;
    00625                         for(int z=1;z<=nrExp;z++)
    00626                         {
    00627                                 pLsig_real=pLsig_real+(ex[z].par[index]-pLmean[j])*(ex[z].par[index]-pLmean[j]);
    00628                         }
    00629                         pLsig_real=sqrt(pLsig_real/nrExp);
    00630                         //*dbg<<"ex["<<i<<"].par["<<index<<"]="<<ex[i].par[index]<<" pLmean["<<j<<"]="<<pLmean[j]<<" pLsig["<<j<<"]="<<pLsig[j]<<" pLsig_real="<<pLsig_real<<endl;
    00631                         // Ma
    00632                         for(int l=1;l<=nrExp;l++)
    00633                         {
    00634                                 Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=1/(pLsig[j]*nrExp);
    00635                                 if(i == l)
    00636                                 {
    00637                                         Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]-1/pLsig[j];
    00638                                 }
    00639                         }
    00640                         k++;
    00641                 }
    00642         }
    00643         // Me
    00644         for(i=eDimold+1;i<=eDim;i++)
    00645         {
    00646                 for(j=nvarFit+1;j<=fitDim;j++)
    00647                 {
    00648                         if(j==fitDimold+i)
    00649                         {
    00650                                 Me[i][j]=Me[i][j]+1;
    00651                         }
    00652                         if(j==nvarFit+i)
    00653                         {
    00654                                 Me[i][j]=Me[i][j]-1;
    00655                         }
    00656                 }
    00657         }
    00658   }
    00659   // Ende Local_Parameter_Constraints
    00660   //****************************************************************************
    00661
    00662   //regularising
    00663   if(globs->reg==TRUE)
    00664     {
    00665       double *diag,**V,**MMa;
    00666       double cond,maxx;
    00667       double reg=0.;
    00668
    00669       V=dmatrix(1,fitDim,1,fitDim);
    00670       MMa=dmatrix(1,aDim,1,fitDim);
    00671       diag=dvector(1,fitDim);
    00672
    00673        for(i=1;i<=aDim;i++)
    00674         {
    00675           for(j=1;j<=fitDim;j++)
    00676             {
    00677               MMa[i][j]=Ma[i][j];
    00678             }
    00679         }
    00680
    00681       //SVD
    00682       svdcmp(MMa,aDim,fitDim,diag,V);
    00683
    00684       //determine condition number
    00685       *dbg << "\ncondition number : ";
    00686       maxx=diag[1];
    00687       cond=diag[1];
    00688       for(i=1;i<=fitDim;i++)
    00689         {
    00690           if(diag[i]>maxx)
    00691             maxx=diag[i];
    00692           if(diag[i]<cond)
    00693             cond=diag[i];
    00694         }
    00695
    00696       //Is the matrix approximately singular?
    00697       for(i=1;i<=fitDim;i++)
    00698         {
    00699           if(fabs(diag[i]/maxx)<=globs->epsilon)
    00700             {
    00701               diag[i]=0.;
    00702               reg=globs->lambda*maxx;
    00703             }
    00704         }
    00705
    00706       cond*=1./maxx;
    00707       *dbg << cond << endl;
    00708       globs->cond=cond;
    00709
    00710       //writing regularised system back
    00711       for(i=1;i<=aDim;i++)
    00712         {
    00713           for(j=1;j<=fitDim;j++)
    00714             {
    00715               Ma[i][j]=0.;
    00716               for(k=1;k<=fitDim;k++)
    00717                 {
    00718                   Ma[i][j]+=MMa[i][k]*diag[k]*V[j][k];
    00719                 }
    00720               if(i <=fitDim)
    00721                 {
    00722                   if(diag[i]==0.)
    00723                     Ma[i][j]+=reg*V[j][i];
    00724                 }
    00725             }
    00726         }
    00727       if(reg!=0.)
    00728         {
    00729           *dbg << "regularisation active, linear dependency of parameters:\n";
    00730           for(i=1;i<=fitDim;i++)
    00731               {
    00732                 if(diag[i]==0.)
    00733                   *dbg << "line " << i << ": ";
    00734                 for(j=1;j<=fitDim;j++)
    00735                   {
    00736                     if(diag[i]==0.)
    00737                       *dbg << V[j][i] << " ";
    00738                   }
    00739                 if(diag[i]==0.)
    00740                   *dbg << endl;
    00741               }
    00742         }
    00743
    00744       free_dmatrix(MMa,1,aDim,1,fitDim);
    00745       free_dmatrix(V,1,fitDim,1,fitDim);
    00746       free_dvector(diag,1,fitDim);
    00747     }//end if(globs->reg==TRUE)
    00748
    00749 #ifdef PRINTBIGSYS
    00750   *dbg << "\nLeast-squares :\n";
    00751    for(i=1;i<=aDim;i++)
    00752     {
    00753       for(j=1;j<=fitDim;j++)
    00754         {
    00755           *dbg << Ma[i][j] << " ";
    00756         }
    00757       *dbg << "| " << Ra[i] << endl;
    00758     }
    00759    *dbg << "\nEquality constraints :\n";
    00760    for(i=1;i<=eDim;i++)
    00761     {
    00762       for(j=1;j<=fitDim;j++)
    00763         {
    00764           *dbg << Me[i][j] << " ";
    00765         }
    00766       *dbg << "| " << Re[i] << endl;
    00767     }
    00768    *dbg << "\nInquality constraints :\n";
    00769    for(i=1;i<=gDim;i++)
    00770     {
    00771       for(j=1;j<=fitDim;j++)
    00772         {
    00773           *dbg << Mg[i][j] << " ";
    00774         }
    00775       *dbg << "| " << Rg[i] << endl;
    00776     }
    00777    dbg->flush();
    00778   #endif
    00779   //allocate solution vector
    00780   dX=(double*)malloc((fitDim+10)*sizeof(double));
    00781   //allocate covariance matrix
    00782   globs->fitdim=fitDim;
    00783   if(globs->covar==NULL)
    00784     globs->covar=dmatrix(1,fitDim,1,fitDim);
    00785
    00786   if(globs->minimiser==1)
    00787     {
    00788       // ###############
    00789       // #     LSEI    #
    00790       // ###############
    00791
    00792       //allocation/initialisation for LSEI
    00793       double *w;
    00794       double *ws;
    00795       double *prgopt;
    00796       double rnorme,rnorml;
    00797       long mdw=aDim+eDim+gDim;
    00798       long mode;
    00799       long ipstor=mdw+npar+2*fitDim+2;
    00800       long wsstor=3*(mdw+npar)+(mdw+npar+2)*(fitDim+8)+fitDim+npar;
    00801       long *ip;
    00802
    00803       w=(double*)malloc(((fitDim+1)*mdw)*sizeof(double));
    00804       ws=(double*)malloc((wsstor+1)*sizeof(double));
    00805       ip=(long*)malloc((ipstor+10)*sizeof(long));
    00806       prgopt=(double*)malloc(30*sizeof(double));
    00807
    00808       xind=0;
    00809       yind=0;
    00810
    00811       for(i=1;i<=fitDim;i++)
    00812         {
    00813           for(j=1;j<=eDim;j++)
    00814             {
    00815               w[xind*mdw+yind]=Me[j][i];
    00816               yind++;
    00817             }
    00818           for(j=1;j<=aDim;j++)
    00819             {
    00820               w[xind*mdw+yind]=Ma[j][i];
    00821               yind++;
    00822             }
    00823           for(j=1;j<=gDim;j++)
    00824             {
    00825               w[xind*mdw+yind]=Mg[j][i];
    00826               yind++;
    00827             }
    00828           yind=0;
    00829           xind++;
    00830         }
    00831
    00832       for(j=1;j<=eDim;j++)
    00833         {
    00834           w[(fitDim)*mdw+yind]=Re[j];
    00835           yind++;
    00836         }
    00837       for(j=1;j<=aDim;j++)
    00838         {
    00839           w[(fitDim)*mdw+yind]=Ra[j];
    00840           yind++;
    00841         }
    00842       for(j=1;j<=gDim;j++)
    00843         {
    00844           w[(fitDim)*mdw+yind]=Rg[j];
    00845           yind++;
    00846         }
    00847
    00848       ip[0]=wsstor;
    00849       ip[1]=ipstor;
    00850
    00851       if(computeCovar==TRUE)
    00852         {
    00853           prgopt[0]=4.0;
    00854           prgopt[1]=1.0;
    00855           prgopt[2]=1.0;
    00856           prgopt[3]=7.0;
    00857           prgopt[4]=10.0;
    00858           prgopt[5]=1.0;
    00859         }
    00860       else
    00861         {
    00862           prgopt[0]=4.0;
    00863           prgopt[1]=0.0;
    00864           prgopt[2]=1.0;
    00865           prgopt[3]=7.0;
    00866           prgopt[4]=10.0;
    00867           prgopt[5]=0.0;
    00868         }
    00869       prgopt[6]=1.0;
    00870
    00871       //call LSEI
    00872       lsei_(w,mdw,eDim,aDim,gDim,fitDim,prgopt,dX,&rnorme,&rnorml,&mode,ws,ip);
    00873
    00874       if(mode>=4)
    00875         {
    00876           cerr << "Minimisation (LSEI) failed.\n";
    00877           throw 1;
    00878         }
    00879
    00880       if(computeCovar==TRUE)
    00881         {
    00882           for (i=1; i<=fitDim; ++i)
    00883             {
    00884               for (j=1;j<=fitDim;++j)
    00885                 globs->covar[i][j] = w[(j-1)*mdw+i-1];
    00886             }
    00887         }
    00888
    00889       free(w);
    00890       free(ws);
    00891       free(ip);
    00892       free(prgopt);
    00893
    00894       // #######  END OF LSEI ######
    00895     }
    00896   else if(globs->minimiser==2)
    00897     {
    00898       //to solve ----->
    00899       //-------------------------
    00900       //
    00901       //  Ma*dX  ~  Ra
    00902       //  Me*dX  =  Re
    00903       //  Mg*dX >=  Rg
    00904       //
    00905       //-------------------------
    00906
    00907
    00908       //initialisation of  E04NCF
    00909       long NCLIN=eDim+gDim;
    00910       double *C=(double*)malloc(((eDim+gDim+1)*fitDim)*sizeof(double));
    00911       double *BL=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double));
    00912       double *BU=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double));
    00913       double *CVEC=(double*)malloc((fitDim+1)*sizeof(double));
    00914       long *ISTATE=(long*)malloc((fitDim+eDim+gDim+1)*sizeof(long));
    00915       long *KX=(long*)malloc((fitDim+1)*sizeof(long));
    00916       double *X=(double*)malloc((1+fitDim)*sizeof(double));
    00917       double *A=(double*)malloc((1+fitDim*aDim)*sizeof(double));
    00918       double *B=(double*)malloc((1+aDim)*sizeof(double));
    00919       long ITER;
    00920       double OBJ;
    00921       double *CLAMDA=(double*)malloc((1+fitDim+eDim+gDim)*sizeof(double));
    00922       long *IWORK=(long*)malloc((1+fitDim)*sizeof(long));
    00923       long LWORK=2*fitDim*fitDim+9*fitDim+6*(eDim+gDim);
    00924       double *WORK=(double*)malloc((1+LWORK)*sizeof(double));
    00925       long IFAIL=-1;
    00926       double **CTMP=dmatrix(1,1+eDim+gDim,1,1+fitDim);
    00927
    00928       //right-hand side
    00929       for(i=1;i<=aDim;i++)
    00930         B[i-1]=Ra[i];
    00931
    00932       //bounds
    00933       for(i=1;i<=fitDim;i++)
    00934         {
    00935           BL[i-1]=-1e25;
    00936           BU[i-1]=1e25;
    00937           X[i-1]=0.;  //initialisation of X
    00938         }
    00939       for(i=fitDim+1;i<=fitDim+eDim;i++)
    00940         {
    00941           BL[i-1]=Re[i-fitDim];
    00942           BU[i-1]=Re[i-fitDim];
    00943         }
    00944       for(i=fitDim+eDim+1;i<=fitDim+eDim+gDim;i++)
    00945         {
    00946           BL[i-1]=Rg[i-fitDim-eDim];
    00947           BU[i-1]=1e25;
    00948         }
    00949
    00950       for(i=1;i<=fitDim;i++)
    00951         {
    00952           for(j=1;j<=eDim+gDim;j++)
    00953             {
    00954               if(j<=eDim)
    00955                 CTMP[j][i]=Me[j][i];
    00956               else
    00957                 CTMP[j][i]=Mg[j-eDim][i];
    00958             }
    00959         }
    00960
    00961       for(i=1;i<=fitDim+eDim+gDim;i++)
    00962       {
    00963         ISTATE[i]=0;
    00964       }
    00965
    00966       cmat2fvec(Ma,aDim,fitDim,A);
    00967       cmat2fvec(CTMP,eDim+gDim,fitDim,C);
    00968
    00969       setopte04ncf_();
    00970       e04ncf_(aDim,fitDim,NCLIN,NCLIN,aDim,C,BL,BU,CVEC,ISTATE,
    00971               KX,X,A,B,ITER,OBJ,CLAMDA,IWORK,fitDim,WORK,LWORK,IFAIL);
    00972       for(i=0;i< fitDim;i++)
    00973       {
    00974         dX[i]=X[i];
    00975       }
    00976
    00977       //free memory
    00978       free(C);
    00979       free(BL);
    00980       free(BU);
    00981       free(CVEC);
    00982       free(ISTATE);
    00983       free(KX);
    00984       free(X);
    00985       free(A);
    00986       free(B);
    00987       free(CLAMDA);
    00988       free(IWORK);
    00989       free(WORK);
    00990       free_dmatrix(CTMP,1,eDim+gDim,1,fitDim);
    00991     }
    00992
    00993   //decondense experiments
    00994   decondense(globs,ex,dX);
    00995
    00996 #ifdef PRINTSTEP
    00997   for(i=1;i<=nrExp;i++)
    00998     {
    00999       *dbg << "\nExperiment #" << i << ":\n\n";
    01000       for(k=1;k<=ex[i].nPoints-1;k++)
    01001         {
    01002           *dbg << "dS(msP=" << k << ") = ";
    01003           for(j=1;j<=nvar;j++)
    01004             *dbg << ex[i].dS[k][j] << " ";
    01005           *dbg << endl;
    01006         }
    01007       *dbg << "dP = ";
    01008       for(j=1;j<=npar;j++)
    01009         *dbg << ex[i].dP[j] << " ";
    01010     }
    01011   *dbg << endl;
    01012   dbg->flush();
    01013 #endif
    01014
    01015
    01016   //free memory
    01017   free_dmatrix(Ma,0,aDim,1,fitDim);
    01018   free_dmatrix(Me,0,eDim,1,fitDim);
    01019   free_dmatrix(Mg,0,gDim,1,fitDim);
    01020   free_dvector(Ra,0,aDim);
    01021   free_dvector(Re,0,eDim);
    01022   free_dvector(Rg,0,gDim);
    01023
    01024   dX[0];
    01025   free(dX);
    01026 }

void xerr_ | ( | char * | msg,
|  | int & | nerr,
|  | int & | iert,
|  | int & | ni,
|  | int & | i1,
|  | int & | i2,
|  | int & | nr,
|  | double & | r1,
|  | double & | r2
| ) |

| Definition at line 1050 of file solvLin.cc.
    01052 {
    01053   cerr << msg << endl;
    01054   if (ni > 0)
    01055     cerr << i1 << " ";
    01056   if (ni > 1)
    01057     cerr << i2 << " ";
    01058   if (nr > 0)
    01059     cerr << r1 << " ";
    01060   if (nr > 1)
    01061     cerr << r2 << " ";
    01062   if (ni > 0 || nr > 0)
    01063     cerr << endl;
    01064   if (iert >= 2)
    01065     cerr << endl;
    01066 }