\section{/home/peifer/diffit/modules/solv\-Lin.cc File Reference}
\label{solvLin_8cc}\index{/home/peifer/diffit/modules/solvLin.cc@{/home/peifer/diffit/modules/solvLin.cc}}
{\tt \#include $<$iostream$>$}\par
{\tt \#include $<$fstream$>$}\par
{\tt \#include $<$math.h$>$}\par
{\tt \#include $<$stdlib.h$>$}\par
{\tt \#include \char`\"{}../def.h\char`\"{}}\par
{\tt \#include \char`\"{}../model.h\char`\"{}}\par
{\tt \#include \char`\"{}../nr.h\char`\"{}}\par
\subsection*{Defines}
\begin{CompactItemize}
\item 
\#define {\bf PRINTBIGSYS}
\end{CompactItemize}
\subsection*{Functions}
\begin{CompactItemize}
\item 
int {\bf lsei\_\-} (double $\ast$w, long \&mdw, long \&me, long \&ma, long \&mg, long \&n, double $\ast$prgopt, double $\ast$x, double $\ast$rnorme, double $\ast$rnorml, long $\ast$mode, double $\ast$ws, long $\ast$ip)
\item 
int {\bf e04ncf\_\-} (long \&M, long \&N, long \&NCLIN, long \&LDC, long \&LDA, double $\ast$C, double $\ast$BL, double $\ast$BU, double $\ast$CVEC, long $\ast$ISTATE, long $\ast$KX, double $\ast$X, double $\ast$A, double $\ast$B, long \&ITER, double \&OBJ, double $\ast$CLAMBDA, long $\ast$IWORK, long \&LIWORK, double $\ast$WORK, long \&LWORK, long \&IFAIL)
\item 
int {\bf setopte04ncf\_\-} ()
\item 
void {\bf condense} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex)
\begin{CompactList}\small\item\em Condenses each experiment. \item\end{CompactList}\item 
void {\bf decondense} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex, double $\ast$d\-X)
\begin{CompactList}\small\item\em Decondenses each experiment. \item\end{CompactList}\item 
void {\bf solv\-Lin} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex, int compute\-Covar)
\item 
void {\bf myerr\_\-} (int $\ast$messg, int \&nmessg)
\item 
void {\bf xerr\_\-} (char $\ast$msg, int \&nerr, int \&iert, int \&ni, int \&i1, int \&i2, int \&nr, double \&r1, double \&r2)
\end{CompactItemize}


\subsection{Define Documentation}
\index{solvLin.cc@{solv\-Lin.cc}!PRINTBIGSYS@{PRINTBIGSYS}}
\index{PRINTBIGSYS@{PRINTBIGSYS}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}\#define PRINTBIGSYS}\label{solvLin_8cc_825af7cb8b39e590564ec13b26ae87af}




Definition at line 14 of file solv\-Lin.cc.

\subsection{Function Documentation}
\index{solvLin.cc@{solv\-Lin.cc}!condense@{condense}}
\index{condense@{condense}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void condense ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex})}\label{solvLin_8cc_988c77101e2bdbd16fecb4084ad4dae3}


Condenses each experiment. 



Definition at line 40 of file solv\-Lin.cc.

References dmatrix(), Glob::elastic, Glob\-Exp::me, Glob\-Exp::mg, Glob\-Exp::nobs, Glob::npar, Glob\-Exp::n\-Points, Glob::nr\-Exp, and Glob\-Exp::nvar.

Referenced by solv\-Lin().

\footnotesize\begin{verbatim}42 {
43   long nExp,msP,i,j,k,l;
44   long nMeas,nvar,nP;
45   long me,mg,nobs,nms;
46   
47   //weakens continuty constraints
48   double elastic=globs->elastic;
49   double **Ea,**Ee,**Eg; //save the previous E-matrices
50 
51   nP=globs->npar;
52   for(nExp=1;nExp<=globs->nrExp;++nExp) 
53     { 
54       nMeas=ex[nExp].nMeasure;
55       nvar=ex[nExp].nvar;
56       me=ex[nExp].me;
57       mg=ex[nExp].mg;
58       nobs=ex[nExp].nobs;
59       nms=ex[nExp].nPoints;
60       nms--;
61       Ea=dmatrix(1,nMeas*nobs,1,nvar);
62       Ee=dmatrix(1,me,1,nvar);
63       Eg=dmatrix(1,mg,1,nvar);
64 
65       //initialisation of backward recursion
66       //------------------------------------
67       //least squares:
68       for(i=1;i<=nMeas;i++)
69         {
70           for(j=1;j<=nobs;j++)
71             {
72               //ua 
73               ex[nExp].ua[(j-1)*nMeas+i]=ex[nExp].residues[i][j]/ex[nExp].sigma[i][j];
74               //Ea
75               for(k=1;k<=nvar;k++)
76                 {
77                   //find out which datapoints are falling in the last MS-interval
78                   //the other derivatives are zero
79                   if(ex[nExp].xMeasure[i] >= ex[nExp].mesh[nms] && 
80                      ex[nExp].xMeasure[i] <= ex[nExp].mesh[nms+1])
81                     {
82                       ex[nExp].Ea[(j-1)*nMeas+i][k]=ex[nExp].dmds[i][j][k]/ex[nExp].sigma[i][j];
83                     }
84                   else
85                     {
86                       ex[nExp].Ea[(j-1)*nMeas+i][k]=0.0;
87                     }
88                   Ea[(j-1)*nMeas+i][k]=ex[nExp].Ea[(j-1)*nMeas+i][k];
89                 }
90               //Pa
91               for(k=1;k<=nP;k++)
92                 {
93                   ex[nExp].Pa[(j-1)*nMeas+i][k]=ex[nExp].dmdp[i][j][k]/ex[nExp].sigma[i][j];
94                 }
95             }
96         }// end of least squares init.
97       
98       //equality constraints
99       for(i=1;i<=me;i++)
100         {
101           //ue
102           ex[nExp].ue[i]=ex[nExp].r2[i];
103           //Ee
104           for(j=1;j<=nvar;j++)
105             {
106               ex[nExp].Ee[i][j]=ex[nExp].dR2ds[i][nms][j];
107               Ee[i][j]=ex[nExp].Ee[i][j];
108             }
109           //Pe
110           for(j=1;j<=nP;j++)
111             ex[nExp].Pe[i][j]=ex[nExp].dR2dp[i][j];
112         }//end of equality constraints
113 
114       //inequality constraints
115       for(i=1;i<=mg;i++)
116         {
117           //ug
118           ex[nExp].ug[i]=ex[nExp].r3[i];
119           //Ee
120           for(j=1;j<=nvar;j++)
121             {
122               ex[nExp].Eg[i][j]=ex[nExp].dR3ds[i][nms][j];
123               Eg[i][j]=ex[nExp].Eg[i][j];
124             }
125           //Pe
126           for(j=1;j<=nP;j++)
127             ex[nExp].Pg[i][j]=ex[nExp].dR3dp[i][j];
128         }//end of inequality constraints
129       //initialisation of backward recursion complete
130     
131       
132       //backward recursion 
133       //-----------------
134       
135       //loop over all multiple shooting intervals
136       for(msP=nms;msP>=2;msP--)
137         {
138           //least squares
139           
140           for(i=1;i<=nMeas;i++)
141             {
142               for(j=1;j<=nobs;j++)
143                 {
144                   //ua 
145                   for(k=1;k<=nvar;k++)
146                     ex[nExp].ua[(j-1)*nMeas+i]+=ex[nExp].Ea[(j-1)*nMeas+i][k]*ex[nExp].h[msP][k]*elastic;
147                   //Ea
148                   for(k=1;k<=nvar;k++)
149                     {
150                       //find out which datapoints are falling in the MS-interval
151                       //the other derivatives are zero
152                       if(ex[nExp].xMeasure[i] >= ex[nExp].mesh[msP-1] && 
153                          ex[nExp].xMeasure[i] <= ex[nExp].mesh[msP])
154                         {
155                           ex[nExp].Ea[(j-1)*nMeas+i][k]=ex[nExp].dmds[i][j][k]/ex[nExp].sigma[i][j];
156                         }
157                       else
158                         {
159                           ex[nExp].Ea[(j-1)*nMeas+i][k]=0.0;
160                         }
161                       for(l=1;l<=nvar;l++)
162                         ex[nExp].Ea[(j-1)*nMeas+i][k]+=Ea[(j-1)*nMeas+i][l]*ex[nExp].dyds[msP][l][k];
163                     }
164                   //Pa
165                   for(k=1;k<=nP;k++)
166                     {
167                       for(l=1;l<=nvar;l++)
168                         {
169                           ex[nExp].Pa[(j-1)*nMeas+i][k]+=Ea[(j-1)*nMeas+i][l]*ex[nExp].dydp[msP][l][k];
170                         }
171                     }
172                 }
173             } //end of least squares
174 
175           //equality constraints
176           for(i=1;i<=me;i++)
177             {
178               //ue
179               for(k=1;k<=nvar;k++)
180                 ex[nExp].ue[i]+=Ee[i][k]*ex[nExp].h[msP][k]*elastic;
181               //Ee
182               for(j=1;j<=nvar;j++)
183                 {
184                   ex[nExp].Ee[i][j]=ex[nExp].dR2ds[i][msP-1][j];
185                     for(l=1;l<=nvar;l++)
186                       ex[nExp].Ee[i][j]+=Ee[i][l]*ex[nExp].dyds[msP][l][j];
187                 }
188               //Pe
189               for(j=1;j<=nP;j++)
190                 {
191                   for(l=1;l<=nvar;l++)
192                     ex[nExp].Pe[i][j]+=Ee[i][l]*ex[nExp].dydp[msP][l][j];
193                 }
194             }//end of equality constraints
195           
196           //inequality constraints
197           for(i=1;i<=mg;i++)
198             {
199               //ue
200               for(k=1;k<=nvar;k++)
201                 ex[nExp].ug[i]+=Eg[i][k]*ex[nExp].h[msP][k]*elastic;
202               //Ee
203               for(j=1;j<=nvar;j++)
204                 {
205                   ex[nExp].Eg[i][j]=ex[nExp].dR3ds[i][msP-1][j];
206                     for(l=1;l<=nvar;l++)
207                       ex[nExp].Eg[i][j]+=Eg[i][l]*ex[nExp].dyds[msP][l][j];
208                 }
209               //Pe
210               for(j=1;j<=nP;j++)
211                 {
212                   for(l=1;l<=nvar;l++)
213                     ex[nExp].Pg[i][j]+=Eg[i][l]*ex[nExp].dydp[msP][l][j];
214                 }
215             }//end of inequality constraints
216 
217           //updating new Ea,Ee,Ea 
218           for(i=1;i<=nvar;i++)
219             {
220               for(j=1;j<=nMeas*nobs;j++)
221                 Ea[j][i]=ex[nExp].Ea[j][i];
222               for(j=1;j<=me;j++)
223                 Ee[j][i]=ex[nExp].Ee[j][i];
224               for(j=1;j<=mg;j++)
225                 Eg[j][i]=ex[nExp].Eg[j][i];           
226             }
227 
228         }//end loop over multiple shooting intervals
229 
230       //clearing memory
231       free_dmatrix(Ea,1,nMeas*nobs,1,nvar);
232       free_dmatrix(Ee,1,me,1,nvar);
233       free_dmatrix(Eg,1,mg,1,nvar);
234     }//end loop over all experiments
235 }
\end{verbatim}\normalsize 


\index{solvLin.cc@{solv\-Lin.cc}!decondense@{decondense}}
\index{decondense@{decondense}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void decondense ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex}, double $\ast$ {\em d\-X})}\label{solvLin_8cc_a5d4e0f36204950ee69a26d51943618c}


Decondenses each experiment. 



Definition at line 238 of file solv\-Lin.cc.

References Glob\-Exp::d\-S, Glob::elastic, FALSE, Glob::npar, Glob\-Exp::n\-Points, Glob::nr\-Exp, Glob\-Exp::nvar, and Glob::y0fix.

\footnotesize\begin{verbatim}240 {
241   long nExp,msP,i,j,k,l;
242   long nvar,npar;
243   long nms,ind=0;
244 
245   //weakens continuty constraints
246   double elastic=globs->elastic;
247 
248   npar=globs->npar;     
249   nvar=ex[1].nvar;
250   
251   //initialise
252   for(nExp=1;nExp<=globs->nrExp;++nExp) 
253     {
254       for(i=1;i<=ex[nExp].nPoints;i++)
255         for(j=1;j<=nvar;j++)
256           ex[nExp].dS[i][j]=0.;
257       for(i=1;i<=npar;i++)
258         ex[nExp].dP[i]=0.;
259 
260       for(i=1;i<=nvar;i++)
261         {
262           if(globs->y0fix[i]!=FALSE)
263             {
264               ex[nExp].dS[1][i]=dX[ind];
265               ind++;
266             }
267         }
268       for(i=1;i<=npar;i++)
269         {
270           if(globs->doP[i]=='L')
271             {
272               ex[nExp].dP[i]=dX[ind];
273               ind++;
274             }
275         }
276     }
277 
278   for(i=1;i<=npar;i++)
279     {
280       if(globs->doP[i]==TRUE)
281         {
282           for(nExp=1;nExp<=globs->nrExp;++nExp) 
283             {
284               ex[nExp].dP[i]=dX[ind];
285             }
286           ind++;
287         }
288     } 
289   
290   //end initialise
291   
292   //forward iteration
293   for(nExp=1;nExp<=globs->nrExp;++nExp) 
294     { 
295       nms=ex[nExp].nPoints;
296       nms--;
297       
298       for(msP=2;msP<=nms;msP++)
299         {
300           for(i=1;i<=nvar;i++)
301             {
302               ex[nExp].dS[msP][i]=elastic*ex[nExp].h[msP][i];
303               for(j=1;j<=nvar;j++)
304                 ex[nExp].dS[msP][i]+=ex[nExp].dyds[msP][i][j]*ex[nExp].dS[msP-1][j];
305               for(j=1;j<=npar;j++)
306                 ex[nExp].dS[msP][i]+=ex[nExp].dydp[msP][i][j]*ex[nExp].dP[j];
307             }
308         }
309       
310     }//end loop over experiments
311 }
\end{verbatim}\normalsize 


\index{solvLin.cc@{solv\-Lin.cc}!e04ncf_@{e04ncf\_\-}}
\index{e04ncf_@{e04ncf\_\-}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}int e04ncf\_\- (long \& {\em M}, long \& {\em N}, long \& {\em NCLIN}, long \& {\em LDC}, long \& {\em LDA}, double $\ast$ {\em C}, double $\ast$ {\em BL}, double $\ast$ {\em BU}, double $\ast$ {\em CVEC}, long $\ast$ {\em ISTATE}, long $\ast$ {\em KX}, double $\ast$ {\em X}, double $\ast$ {\em A}, double $\ast$ {\em B}, long \& {\em ITER}, double \& {\em OBJ}, double $\ast$ {\em CLAMBDA}, long $\ast$ {\em IWORK}, long \& {\em LIWORK}, double $\ast$ {\em WORK}, long \& {\em LWORK}, long \& {\em IFAIL})}\label{solvLin_8cc_8683b6d97ba48de409e80c13664d1ebd}


\index{solvLin.cc@{solv\-Lin.cc}!lsei_@{lsei\_\-}}
\index{lsei_@{lsei\_\-}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}int lsei\_\- (double $\ast$ {\em w}, long \& {\em mdw}, long \& {\em me}, long \& {\em ma}, long \& {\em mg}, long \& {\em n}, double $\ast$ {\em prgopt}, double $\ast$ {\em x}, double $\ast$ {\em rnorme}, double $\ast$ {\em rnorml}, long $\ast$ {\em mode}, double $\ast$ {\em ws}, long $\ast$ {\em ip})}\label{solvLin_8cc_1d5f422031db31661e93ce27e6efdd16}


\index{solvLin.cc@{solv\-Lin.cc}!myerr_@{myerr\_\-}}
\index{myerr_@{myerr\_\-}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void myerr\_\- (int $\ast$ {\em messg}, int \& {\em nmessg})}\label{solvLin_8cc_f7817de873b0d1db29c8e0af563e15ec}




Definition at line 1034 of file solv\-Lin.cc.

\footnotesize\begin{verbatim}1035 {
1036   char *s = (char *) messg;
1037   cerr << endl;
1038   for (int i = 0; i < nmessg; i++)
1039     cerr << s[i];
1040   cerr << endl;
1041 }
\end{verbatim}\normalsize 


\index{solvLin.cc@{solv\-Lin.cc}!setopte04ncf_@{setopte04ncf\_\-}}
\index{setopte04ncf_@{setopte04ncf\_\-}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}int setopte04ncf\_\- ()}\label{solvLin_8cc_fda1e06079b0e35830652847fda3ecd7}


\index{solvLin.cc@{solv\-Lin.cc}!solvLin@{solvLin}}
\index{solvLin@{solvLin}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void solv\-Lin ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex}, int {\em compute\-Covar})}\label{solvLin_8cc_ff7b6c2bcde7c9861871a4f9df066a07}




Definition at line 313 of file solv\-Lin.cc.

References condense(), FALSE, Glob::npar, Glob::nr\-Exp, Glob\-Exp::nvar, and Glob::y0fix.

\footnotesize\begin{verbatim}314 {
315   
316   long i,j,k,l,xind=0,yind=0,nL;
317   long nvarFit=0;  // # der zu fittenden Anfangswerte
318   long nLpara=0;   // # lokale Parameter
319   long nparL=0;    // nLpara*nrExp
320   long npar=globs->npar;
321   long nvar=ex[1].nvar;
322   long nrExp=globs->nrExp;
323   long nparG=0;
324   long fitDim=0;
325   long aDim=0,eDim=0,gDim=0;
326 
327   //big system (condensed)
328   double **Ma,*Ra;
329   double **Me,*Re;
330   double **Mg,*Rg;
331 
332   //solution
333   double *dX;
334   //condense experiments
335   condense(globs,ex);
336 
337   //determine dimension of big system
338   for(i=1;i<=nvar;i++)
339     {
340       if(globs->y0fix[i]!=FALSE)
341         nvarFit++;
342     }
343   nvarFit*=nrExp;
344   for(i=1;i<=npar;i++)
345     {
346       if(globs->doP[i]=='L')
347         nparL++;
348       else if(globs->doP[i]!=FALSE)
349         nparG++;
350     }
351   nLpara=nparL;
352   nparL*=nrExp;
353   fitDim=nvarFit+nparL+nparG;
354   for(i=1;i<=nrExp;i++)
355     {
356       aDim+=ex[i].nMeasure*ex[i].nobs;
357       eDim+=ex[i].me;
358       gDim+=ex[i].mg;
359     }
360   if(fitDim==0)
361     {
362       cerr << "Nothing to fit.\n";
363       throw 1;
364     }
365 #ifdef PRINTDIMENSION
366   //print dimension of bigsys;
367   *dbg << "Dimension of Big-System\n";
368   *dbg << "-----------------------\n";
369   *dbg << "Least-squares:\n";
370   *dbg << "#Rows : " << aDim  << "  #Columns : " << fitDim << endl;
371   *dbg << "Equality constraints:\n";
372   *dbg << "#Rows : " << eDim  << "  #Columns : " << fitDim << endl;  
373   *dbg << "Inequality constraints:\n";
374   *dbg << "#Rows : " << gDim  << "  #Columns : " << fitDim;  
375   *dbg << "\n\n";
376   dbg->flush();
377 #endif
378 
379 
380   long aDimold=aDim;
381   long eDimold=eDim;
382   long fitDimold=fitDim;
383 
384   // Beginn Local_Parameter_Constraints
385   if(globs->faktorLexist)
386   {
387         // Erweiterung der Dimensionen um nLpara*nExp
388         aDim=aDim+nparL;
389         eDim=eDim+nparL;
390         fitDim=fitDim+nparL;
391   }
392   // Ende Local_Parameter_Constraints
393   
394   //allocate memory
395   Ma=dmatrix(0,aDim,1,fitDim);
396   Me=dmatrix(0,eDim,1,fitDim);
397   Mg=dmatrix(0,gDim,1,fitDim);
398   Ra=dvector(0,aDim);
399   Re=dvector(0,eDim);
400   Rg=dvector(0,gDim);
401 
402   //initialise 
403   for(i=1;i<=aDim;i++)
404     {
405       Ra[i]=0.;
406       for(j=1;j<=fitDim;j++)
407         Ma[i][j]=0.;
408     }
409    for(i=1;i<=eDim;i++)
410     {
411       Re[i]=0.;
412       for(j=1;j<=fitDim;j++)
413         Me[i][j]=0.;
414     } 
415    for(i=1;i<=gDim;i++)
416      {
417       Rg[i]=0.;
418       for(j=1;j<=fitDim;j++)
419         Mg[i][j]=0.;
420     } 
421 
422   //least-squares
423   for(i=1;i<=nrExp;i++)
424     {
425       for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)    
426         Ra[yind+k]=-ex[i].ua[k];
427       //initial values
428       for(j=1;j<=nvar;j++)
429         {
430           if(globs->y0fix[j]!=FALSE)
431             {
432               xind++;
433               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)    
434                 Ma[yind+k][xind]=ex[i].Ea[k][j];
435             }
436           
437         }
438       //parameters
439       nL=0;
440       for(j=1;j<=npar;j++)
441         {
442           if(globs->doP[j]=='L')
443             {
444               xind++;
445               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
446                 Ma[yind+k][xind]=ex[i].Pa[k][j];
447             }
448           else if(globs->doP[j]!=FALSE)
449             {
450               nL++;
451               for(k=1;k<=ex[i].nMeasure*ex[i].nobs;k++)
452                 Ma[yind+k][nvarFit+nparL+nL]=ex[i].Pa[k][j];
453             }
454         }
455       yind=yind+(ex[i].nMeasure*ex[i].nobs);
456     }//end loop over experiments
457   //-------------> end lest squares
458   
459   //equality constraints
460   yind=0;
461   xind=0;
462   for(i=1;i<=nrExp;i++)
463     {
464       for(k=1;k<=ex[i].me;k++)    
465         Re[yind+k]=-ex[i].ue[k];
466       //initial values
467       for(j=1;j<=nvar;j++)
468         {
469           if(globs->y0fix[j]!=FALSE)
470             {
471               xind++;
472               for(k=1;k<=ex[i].me;k++)    
473                 Me[yind+k][xind]=ex[i].Ee[k][xind];
474             }
475         }
476       //parameters
477       nL=0;
478       for(j=1;j<=npar;j++)
479         {
480           if(globs->doP[j]=='L')
481             {
482               xind++;
483               for(k=1;k<=ex[i].me;k++)
484                 Me[yind+k][xind]=ex[i].Pe[k][j];
485             }
486           else if(globs->doP[j]!=FALSE)
487             {
488               nL++;
489               for(k=1;k<=ex[i].me;k++)
490                 Me[yind+k][nvarFit+nparL+nL]=ex[i].Pe[k][j];
491             }
492         }
493       yind+=ex[i].me;
494     }//end loop over experiments
495   //-------------> end equality  constraints
496 
497   //inequality constraints
498   yind=0;
499   xind=0;
500   for(i=1;i<=nrExp;i++)
501     {
502       for(k=1;k<=ex[i].mg;k++)    
503         Rg[yind+k]=-ex[i].ug[k];
504       //initial values
505       for(j=1;j<=nvar;j++)
506         {
507           if(globs->y0fix[j]!=FALSE)
508             {
509               xind++;
510               for(k=1;k<=ex[i].mg;k++)  
511                 Mg[yind+k][xind]=ex[i].Eg[k][xind];
512             }
513         }
514       //parameters
515       nL=0;
516       for(j=1;j<=npar;j++)
517         {
518           if(globs->doP[j]=='L')
519             {
520               xind++;
521               for(k=1;k<=ex[i].mg;k++)
522                 Mg[yind+k][xind]=ex[i].Pg[k][j];
523             }
524           else if(globs->doP[j]!=FALSE)
525             {
526               nL++;
527               for(k=1;k<=ex[i].mg;k++)
528                 Mg[yind+k][nvarFit+nparL+nL]=ex[i].Pg[k][j];
529             }
530         }
531       yind+=ex[i].mg;
532     }//end loop over experiments
533   //-------------> end equality  constraints
534   
535   //Big System Ready to Solve
536   //-------------------------
537   //
538   //  Ma*dX  ~  Ra
539   //  Me*dX  =  Re
540   //  Mg*dX >=  Rg
541   //
542   //-------------------------
543   
544   //****************************************************************************
545   // Beginn Local_Parameter_Constraints
546   // Erweiterung der Systems auf Local_Parameter_Constraints:
547   //
548   // pLmean und pLsig bestimmen
549   if(globs->faktorLexist)
550   {
551         double *pLmean, *pLsig;
552         int index;
553         int index_local_parameter;
554         pLmean=dvector(1,nLpara);
555         pLsig=dvector(1,nLpara);
556   
557 //      for(i=1;i<=nLpara;i++)
558 //      {
559 //         cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl;
560 //         if(globs->faktorL[i]<0)
561 //         {
562 //            //i;
563 //            break;
564 //         }
565 //      }
566 //      cerr<<"i="<< i<<" faktorL[i]="<<globs->faktorL[i]<<endl;
567 //      if(i<=nLpara)
568 //      {
569 //         cerr<<"number of local parameter constraints too small"<<endl;
570 //         throw(1);
571 //      }
572         for(i=1;i<=nLpara;i++)
573         {
574                 pLmean[i]=0;
575                 index=0;
576                 index_local_parameter=0;
577                 while(index_local_parameter!=i)
578                 {
579                         index++;
580 //                      cerr<<"doP[index]="<<globs->doP[index]<<" index="<<index<<" i="<<i<<endl;
581                         if(globs->doP[index]=='L')
582                         {
583                                 index_local_parameter++;
584 //                              cerr<<"index="<<index<<" i="<<i<<endl;
585                         }
586                         
587                 }               
588                 for(j=1;j<=nrExp;j++)
589                 {
590                         pLmean[i]=pLmean[i]+ex[j].par[index]; 
591 //                      cerr<<"ex[j].par[index]="<<ex[j].par[index]<<endl;
592                 }
593                 pLmean[i]=pLmean[i]/double(nrExp);
594           //    cerr<<"pLmean["<<i<<"]="<<pLmean[i]<<" globs->faktorL["<<i<<"]"<<globs->faktorL[i]<<endl;
595         }
596         // Berechnung pLsig
597         for(i=1;i<=nLpara;i++)
598         {
599                 pLsig[i]=(globs->faktorL[i]*pLmean[i]-pLmean[i])/(1+globs->faktorL[i]);
600                 if(pLsig[i]<0)
601                 {
602                 cerr<<"ERROR - Negative varianz bei den Lokalen Parameter Constraints..."<<endl;
603                 }
604         }       
605         // Ra,Re
606         k=1;
607         for(i=1;i<=nrExp;i++)
608         {
609                 for(j=1;j<=nLpara;j++)
610                 {       
611                         index=0;
612                         index_local_parameter=0;
613                         while(index_local_parameter!=j)
614                         {
615                                 index++;
616                                 if(globs->doP[index]=='L')
617                                 {
618                                         index_local_parameter++;
619                                 }
620                         
621                         }  
622                         Ra[aDimold+k]=(ex[i].par[index]-pLmean[j])/pLsig[j];
623                         Re[eDimold+k]=0;
624                         double pLsig_real=0;
625                         for(int z=1;z<=nrExp;z++)
626                         {
627                                 pLsig_real=pLsig_real+(ex[z].par[index]-pLmean[j])*(ex[z].par[index]-pLmean[j]);
628                         }
629                         pLsig_real=sqrt(pLsig_real/nrExp);
630                         //*dbg<<"ex["<<i<<"].par["<<index<<"]="<<ex[i].par[index]<<" pLmean["<<j<<"]="<<pLmean[j]<<" pLsig["<<j<<"]="<<pLsig[j]<<" pLsig_real="<<pLsig_real<<endl;
631                         // Ma
632                         for(int l=1;l<=nrExp;l++)
633                         {
634                                 Ma[aDimold+(i-1)*nLpara+j][fitDimold+(l-1)*nLpara+j]=1/(pLsig[j]*nrExp);
635                                 if(i == l)
636                                 {
637                                         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];
638                                 }
639                         }
640                         k++;
641                 }
642         }
643         // Me
644         for(i=eDimold+1;i<=eDim;i++)
645         {
646                 for(j=nvarFit+1;j<=fitDim;j++)   
647                 {
648                         if(j==fitDimold+i)
649                         {
650                                 Me[i][j]=Me[i][j]+1;
651                         }
652                         if(j==nvarFit+i)
653                         {
654                                 Me[i][j]=Me[i][j]-1;
655                         }
656                 }
657         }
658   }
659   // Ende Local_Parameter_Constraints
660   //****************************************************************************
661   
662   //regularising
663   if(globs->reg==TRUE)
664     {
665       double *diag,**V,**MMa;
666       double cond,maxx;
667       double reg=0.;
668 
669       V=dmatrix(1,fitDim,1,fitDim);
670       MMa=dmatrix(1,aDim,1,fitDim);
671       diag=dvector(1,fitDim);
672 
673        for(i=1;i<=aDim;i++)
674         {
675           for(j=1;j<=fitDim;j++)
676             {
677               MMa[i][j]=Ma[i][j];
678             }     
679         }
680       
681       //SVD
682       svdcmp(MMa,aDim,fitDim,diag,V);
683       
684       //determine condition number
685       *dbg << "\ncondition number : ";
686       maxx=diag[1];
687       cond=diag[1];
688       for(i=1;i<=fitDim;i++)
689         {
690           if(diag[i]>maxx)
691             maxx=diag[i];
692           if(diag[i]<cond)
693             cond=diag[i];
694         }
695 
696       //Is the matrix approximately singular?
697       for(i=1;i<=fitDim;i++)
698         {
699           if(fabs(diag[i]/maxx)<=globs->epsilon)
700             {
701               diag[i]=0.;
702               reg=globs->lambda*maxx;
703             }
704         }
705       
706       cond*=1./maxx;
707       *dbg << cond << endl;
708       globs->cond=cond;
709       
710       //writing regularised system back
711       for(i=1;i<=aDim;i++)
712         {
713           for(j=1;j<=fitDim;j++)
714             {
715               Ma[i][j]=0.;
716               for(k=1;k<=fitDim;k++)
717                 {
718                   Ma[i][j]+=MMa[i][k]*diag[k]*V[j][k];
719                 }
720               if(i <=fitDim)
721                 {
722                   if(diag[i]==0.)
723                     Ma[i][j]+=reg*V[j][i];
724                 }
725             }
726         }
727       if(reg!=0.)
728         {
729           *dbg << "regularisation active, linear dependency of parameters:\n";
730           for(i=1;i<=fitDim;i++)
731               {
732                 if(diag[i]==0.)
733                   *dbg << "line " << i << ": ";
734                 for(j=1;j<=fitDim;j++)
735                   {
736                     if(diag[i]==0.)
737                       *dbg << V[j][i] << " ";
738                   }
739                 if(diag[i]==0.)
740                   *dbg << endl;
741               }
742         }
743 
744       free_dmatrix(MMa,1,aDim,1,fitDim);
745       free_dmatrix(V,1,fitDim,1,fitDim);
746       free_dvector(diag,1,fitDim);
747     }//end if(globs->reg==TRUE)
748 
749 #ifdef PRINTBIGSYS
750   *dbg << "\nLeast-squares :\n";
751    for(i=1;i<=aDim;i++)
752     {
753       for(j=1;j<=fitDim;j++)
754         {
755           *dbg << Ma[i][j] << " ";
756         }
757       *dbg << "| " << Ra[i] << endl;
758     }  
759    *dbg << "\nEquality constraints :\n";
760    for(i=1;i<=eDim;i++)
761     {
762       for(j=1;j<=fitDim;j++)
763         {
764           *dbg << Me[i][j] << " ";
765         }
766       *dbg << "| " << Re[i] << endl;
767     }
768    *dbg << "\nInquality constraints :\n";
769    for(i=1;i<=gDim;i++)
770     {
771       for(j=1;j<=fitDim;j++)
772         {
773           *dbg << Mg[i][j] << " ";
774         }
775       *dbg << "| " << Rg[i] << endl;
776     }
777    dbg->flush();
778   #endif
779   //allocate solution vector
780   dX=(double*)malloc((fitDim+10)*sizeof(double));
781   //allocate covariance matrix
782   globs->fitdim=fitDim;
783   if(globs->covar==NULL)
784     globs->covar=dmatrix(1,fitDim,1,fitDim);
785 
786   if(globs->minimiser==1)
787     {
788       // ###############
789       // #     LSEI    #
790       // ###############
791       
792       //allocation/initialisation for LSEI
793       double *w;
794       double *ws;
795       double *prgopt;
796       double rnorme,rnorml;
797       long mdw=aDim+eDim+gDim;
798       long mode;
799       long ipstor=mdw+npar+2*fitDim+2;
800       long wsstor=3*(mdw+npar)+(mdw+npar+2)*(fitDim+8)+fitDim+npar;
801       long *ip;
802       
803       w=(double*)malloc(((fitDim+1)*mdw)*sizeof(double));
804       ws=(double*)malloc((wsstor+1)*sizeof(double));
805       ip=(long*)malloc((ipstor+10)*sizeof(long));
806       prgopt=(double*)malloc(30*sizeof(double));
807       
808       xind=0;
809       yind=0;
810       
811       for(i=1;i<=fitDim;i++)
812         {
813           for(j=1;j<=eDim;j++)
814             {
815               w[xind*mdw+yind]=Me[j][i];
816               yind++;
817             }
818           for(j=1;j<=aDim;j++)
819             {
820               w[xind*mdw+yind]=Ma[j][i];
821               yind++;
822             }
823           for(j=1;j<=gDim;j++)
824             {
825               w[xind*mdw+yind]=Mg[j][i];
826               yind++;
827             } 
828           yind=0;
829           xind++;
830         }
831  
832       for(j=1;j<=eDim;j++)
833         {
834           w[(fitDim)*mdw+yind]=Re[j];
835           yind++;
836         }
837       for(j=1;j<=aDim;j++)
838         {
839           w[(fitDim)*mdw+yind]=Ra[j];
840           yind++;
841         }
842       for(j=1;j<=gDim;j++)
843         {
844           w[(fitDim)*mdw+yind]=Rg[j];
845           yind++;
846         }
847       
848       ip[0]=wsstor;
849       ip[1]=ipstor;
850       
851       if(computeCovar==TRUE)
852         {
853           prgopt[0]=4.0; 
854           prgopt[1]=1.0; 
855           prgopt[2]=1.0;
856           prgopt[3]=7.0;
857           prgopt[4]=10.0;
858           prgopt[5]=1.0;
859         }
860       else
861         {
862           prgopt[0]=4.0; 
863           prgopt[1]=0.0; 
864           prgopt[2]=1.0;
865           prgopt[3]=7.0;
866           prgopt[4]=10.0;
867           prgopt[5]=0.0;
868         }
869       prgopt[6]=1.0;
870       
871       //call LSEI
872       lsei_(w,mdw,eDim,aDim,gDim,fitDim,prgopt,dX,&rnorme,&rnorml,&mode,ws,ip);
873       
874       if(mode>=4)
875         {
876           cerr << "Minimisation (LSEI) failed.\n";
877           throw 1;
878         }
879       
880       if(computeCovar==TRUE)
881         {
882           for (i=1; i<=fitDim; ++i)
883             {
884               for (j=1;j<=fitDim;++j)
885                 globs->covar[i][j] = w[(j-1)*mdw+i-1];
886             } 
887         }
888       
889       free(w);
890       free(ws);
891       free(ip);
892       free(prgopt);
893       
894       // #######  END OF LSEI ######
895     }
896   else if(globs->minimiser==2)
897     {
898       //to solve ----->
899       //-------------------------
900       //
901       //  Ma*dX  ~  Ra
902       //  Me*dX  =  Re
903       //  Mg*dX >=  Rg
904       //
905       //-------------------------
906 
907 
908       //initialisation of  E04NCF
909       long NCLIN=eDim+gDim;
910       double *C=(double*)malloc(((eDim+gDim+1)*fitDim)*sizeof(double));
911       double *BL=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double));
912       double *BU=(double*)malloc((fitDim+eDim+gDim+1)*sizeof(double));
913       double *CVEC=(double*)malloc((fitDim+1)*sizeof(double));
914       long *ISTATE=(long*)malloc((fitDim+eDim+gDim+1)*sizeof(long));
915       long *KX=(long*)malloc((fitDim+1)*sizeof(long));
916       double *X=(double*)malloc((1+fitDim)*sizeof(double));
917       double *A=(double*)malloc((1+fitDim*aDim)*sizeof(double));
918       double *B=(double*)malloc((1+aDim)*sizeof(double));
919       long ITER;
920       double OBJ;
921       double *CLAMDA=(double*)malloc((1+fitDim+eDim+gDim)*sizeof(double));
922       long *IWORK=(long*)malloc((1+fitDim)*sizeof(long));
923       long LWORK=2*fitDim*fitDim+9*fitDim+6*(eDim+gDim);
924       double *WORK=(double*)malloc((1+LWORK)*sizeof(double));
925       long IFAIL=-1;
926       double **CTMP=dmatrix(1,1+eDim+gDim,1,1+fitDim);
927 
928       //right-hand side
929       for(i=1;i<=aDim;i++)
930         B[i-1]=Ra[i];
931       
932       //bounds
933       for(i=1;i<=fitDim;i++)
934         {
935           BL[i-1]=-1e25;
936           BU[i-1]=1e25;
937           X[i-1]=0.;  //initialisation of X
938         }
939       for(i=fitDim+1;i<=fitDim+eDim;i++)
940         {
941           BL[i-1]=Re[i-fitDim];
942           BU[i-1]=Re[i-fitDim];
943         }
944       for(i=fitDim+eDim+1;i<=fitDim+eDim+gDim;i++)
945         {
946           BL[i-1]=Rg[i-fitDim-eDim];
947           BU[i-1]=1e25;
948         }
949       
950       for(i=1;i<=fitDim;i++)
951         {
952           for(j=1;j<=eDim+gDim;j++)
953             {
954               if(j<=eDim)
955                 CTMP[j][i]=Me[j][i];
956               else
957                 CTMP[j][i]=Mg[j-eDim][i];
958             }
959         }
960       
961       for(i=1;i<=fitDim+eDim+gDim;i++)
962       {
963         ISTATE[i]=0;
964       }
965       
966       cmat2fvec(Ma,aDim,fitDim,A);
967       cmat2fvec(CTMP,eDim+gDim,fitDim,C);
968       
969       setopte04ncf_();
970       e04ncf_(aDim,fitDim,NCLIN,NCLIN,aDim,C,BL,BU,CVEC,ISTATE,
971               KX,X,A,B,ITER,OBJ,CLAMDA,IWORK,fitDim,WORK,LWORK,IFAIL);
972       for(i=0;i< fitDim;i++)
973       {
974         dX[i]=X[i];
975       }
976       
977       //free memory
978       free(C);
979       free(BL);
980       free(BU);
981       free(CVEC);
982       free(ISTATE);
983       free(KX);
984       free(X);
985       free(A);
986       free(B);
987       free(CLAMDA);
988       free(IWORK);
989       free(WORK);
990       free_dmatrix(CTMP,1,eDim+gDim,1,fitDim);
991     }
992 
993   //decondense experiments
994   decondense(globs,ex,dX);
995 
996 #ifdef PRINTSTEP
997   for(i=1;i<=nrExp;i++)
998     {
999       *dbg << "\nExperiment #" << i << ":\n\n";
1000       for(k=1;k<=ex[i].nPoints-1;k++)
1001         {
1002           *dbg << "dS(msP=" << k << ") = ";
1003           for(j=1;j<=nvar;j++)
1004             *dbg << ex[i].dS[k][j] << " ";
1005           *dbg << endl;
1006         }
1007       *dbg << "dP = ";
1008       for(j=1;j<=npar;j++)
1009         *dbg << ex[i].dP[j] << " ";
1010     }
1011   *dbg << endl;
1012   dbg->flush();
1013 #endif
1014   
1015 
1016   //free memory
1017   free_dmatrix(Ma,0,aDim,1,fitDim);
1018   free_dmatrix(Me,0,eDim,1,fitDim);
1019   free_dmatrix(Mg,0,gDim,1,fitDim);
1020   free_dvector(Ra,0,aDim);
1021   free_dvector(Re,0,eDim);
1022   free_dvector(Rg,0,gDim);
1023 
1024   dX[0];
1025   free(dX);
1026 }
\end{verbatim}\normalsize 


\index{solvLin.cc@{solv\-Lin.cc}!xerr_@{xerr\_\-}}
\index{xerr_@{xerr\_\-}!solvLin.cc@{solv\-Lin.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void xerr\_\- (char $\ast$ {\em msg}, int \& {\em nerr}, int \& {\em iert}, int \& {\em ni}, int \& {\em i1}, int \& {\em i2}, int \& {\em nr}, double \& {\em r1}, double \& {\em r2})}\label{solvLin_8cc_b524a42b83c6fdc3926032718ff435ce}




Definition at line 1050 of file solv\-Lin.cc.

\footnotesize\begin{verbatim}1052 {
1053   cerr << msg << endl;
1054   if (ni > 0)
1055     cerr << i1 << " ";
1056   if (ni > 1)
1057     cerr << i2 << " ";
1058   if (nr > 0)
1059     cerr << r1 << " ";
1060   if (nr > 1)
1061     cerr << r2 << " ";
1062   if (ni > 0 || nr > 0)
1063     cerr << endl;
1064   if (iert >= 2)
1065     cerr << endl;
1066 }
\end{verbatim}\normalsize 


