Cantera  2.0
BasisOptimize.cpp
Go to the documentation of this file.
1 /**
2  * @file BasisOptimize.cpp
3  * Functions which calculation optimized basis of the
4  * stoichiometric coefficient matrix (see /ref equil functions)
5  */
6 #include "cantera/base/ct_defs.h"
7 #include "cantera/base/global.h"
8 #include "cantera/thermo/ThermoPhase.h"
9 #include "cantera/equil/MultiPhase.h"
10 
11 #include <cstring>
12 
13 using namespace Cantera;
14 using namespace std;
15 
16 #ifdef DEBUG_MODE
17 namespace Cantera
18 {
19 int BasisOptimize_print_lvl = 0;
20 }
21 
22 //! Print a string within a given space limit. This routine limits the amount of the string that will be printed to a
23 //! maximum of "space" characters.
24 /*!
25  *
26  * @param str String -> must be null terminated.
27  * @param space space limit for the printing.
28  * @param alignment 0 centered
29  * 1 right aligned
30  * 2 left aligned
31  */
32 static void print_stringTrunc(const char* str, int space, int alignment);
33 #endif
34 
35 
36 //! Finds the location of the maximum component in a double vector INPUT
37 /*!
38  * @param x Vector to search
39  * @param j j <= i < n : i is the range of indices to search in X(*)
40  * @param n Length of the vector
41  *
42  * @return index of the greatest value on X(*) searched
43  */
44 static size_t amax(double* x, size_t j, size_t n);
45 
46 //! Invert an nxn matrix and solve m rhs's
47 /*!
48  *
49  * Solve C X + B = 0;
50  *
51  * This routine uses Gauss elimination and is optimized for the solution
52  * of lots of rhs's.
53  * A crude form of row pivoting is used here.
54  *
55  * @param c C is the matrix to be inverted
56  * @param idem first dimension in the calling routine
57  * idem >= n must be true
58  * @param n number of rows and columns in the matrix
59  * @param b rhs of the matrix problem
60  * @param m number of rhs to be solved for
61  *
62  * c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number
63  * j = column number
64  * b[i+j*idem] = b_i_j = vectors of rhs's: i = row number
65  * j = column number
66  * (each column is a new rhs)
67  *
68  * @return Retuns the value
69  * 1 : Matrix is singular
70  * 0 : solution is OK
71  *
72  * The solution is returned in the matrix b.
73  */
74 static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m);
75 
76 /*
77  * Choose the optimum basis for the calculations. This is done by
78  * choosing the species with the largest mole fraction
79  * not currently a linear combination of the previous components.
80  * Then, calculate the stoichiometric coefficient matrix for that
81  * basis.
82  *
83  * Calculates the identity of the component species in the mechanism.
84  * Rearranges the solution data to put the component data at the
85  * front of the species list.
86  *
87  * Then, calculates SC(J,I) the formation reactions for all noncomponent
88  * species in the mechanism.
89  *
90  * Input
91  * ---------
92  * mphase Pointer to the multiphase object. Contains the
93  * species mole fractions, which are used to pick the
94  * current optimal species component basis.
95  * orderVectorElement
96  * Order vector for the elements. The element rows
97  * in the formula matrix are
98  * rearranged according to this vector.
99  * orderVectorSpecies
100  * Order vector for the species. The species are
101  * rearranged according to this formula. The first
102  * nCompoments of this vector contain the calculated
103  * species components on exit.
104  * doFormRxn If true, the routine calculates the formation
105  * reaction matrix based on the calculated
106  * component species. If false, this step is skipped.
107  *
108  * Output
109  * ---------
110  * usedZeroedSpecies = If true, then a species with a zero concentration
111  * was used as a component. The problem may be
112  * converged.
113  * formRxnMatrix
114  *
115  * Return
116  * --------------
117  * returns the number of components.
118  *
119  *
120  */
121 size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
122  MultiPhase* mphase, std::vector<size_t>& orderVectorSpecies,
123  std::vector<size_t>& orderVectorElements,
124  vector_fp& formRxnMatrix)
125 {
126 
127  size_t j, jj, k=0, kk, l, i, jl, ml;
128  bool lindep;
129  std::string ename;
130  std::string sname;
131  /*
132  * Get the total number of elements defined in the multiphase object
133  */
134  size_t ne = mphase->nElements();
135  /*
136  * Get the total number of species in the multiphase object
137  */
138  size_t nspecies = mphase->nSpecies();
139  doublereal tmp;
140  doublereal const USEDBEFORE = -1;
141 
142  /*
143  * Perhaps, initialize the element ordering
144  */
145  if (orderVectorElements.size() < ne) {
146  orderVectorElements.resize(ne);
147  for (j = 0; j < ne; j++) {
148  orderVectorElements[j] = j;
149  }
150  }
151 
152  /*
153  * Perhaps, initialize the species ordering
154  */
155  if (orderVectorSpecies.size() != nspecies) {
156  orderVectorSpecies.resize(nspecies);
157  for (k = 0; k < nspecies; k++) {
158  orderVectorSpecies[k] = k;
159  }
160  }
161 
162 #ifdef DEBUG_MODE
163  double molSave = 0.0;
164  if (BasisOptimize_print_lvl >= 1) {
165  writelog(" ");
166  for (i=0; i<77; i++) {
167  writelog("-");
168  }
169  writelog("\n");
170  writelog(" --- Subroutine BASOPT called to ");
171  writelog("calculate the number of components and ");
172  writelog("evaluate the formation matrix\n");
173  if (BasisOptimize_print_lvl > 0) {
174  writelog(" ---\n");
175 
176  writelog(" --- Formula Matrix used in BASOPT calculation\n");
177  writelog(" --- Species | Order | ");
178  for (j = 0; j < ne; j++) {
179  jj = orderVectorElements[j];
180  writelog(" ");
181  ename = mphase->elementName(jj);
182  print_stringTrunc(ename.c_str(), 4, 1);
183  writelogf("(%1d)", j);
184  }
185  writelog("\n");
186  for (k = 0; k < nspecies; k++) {
187  kk = orderVectorSpecies[k];
188  writelog(" --- ");
189  sname = mphase->speciesName(kk);
190  print_stringTrunc(sname.c_str(), 11, 1);
191  writelogf(" | %4d |", k);
192  for (j = 0; j < ne; j++) {
193  jj = orderVectorElements[j];
194  double num = mphase->nAtoms(kk,jj);
195  writelogf("%6.1g ", num);
196  }
197  writelog("\n");
198  }
199  writelog(" --- \n");
200  }
201  }
202 #endif
203 
204  /*
205  * Calculate the maximum value of the number of components possible
206  * It's equal to the minimum of the number of elements and the
207  * number of total species.
208  */
209  size_t nComponents = std::min(ne, nspecies);
210  size_t nNonComponents = nspecies - nComponents;
211  /*
212  * Set this return variable to false
213  */
214  *usedZeroedSpecies = false;
215 
216  /*
217  * Create an array of mole numbers
218  */
219  vector_fp molNum(nspecies,0.0);
220  mphase->getMoles(DATA_PTR(molNum));
221 
222  /*
223  * Other workspace
224  */
225  vector_fp sm(ne*ne, 0.0);
226  vector_fp ss(ne, 0.0);
227  vector_fp sa(ne, 0.0);
228  if (formRxnMatrix.size() < nspecies*ne) {
229  formRxnMatrix.resize(nspecies*ne, 0.0);
230  }
231 
232 #ifdef DEBUG_MODE
233  /*
234  * For debugging purposes keep an unmodified copy of the array.
235  */
236  vector_fp molNumBase(molNum);
237 #endif
238 
239 
240  size_t jr = npos;
241  /*
242  * Top of a loop of some sort based on the index JR. JR is the
243  * current number of component species found.
244  */
245  do {
246  ++jr;
247  /* - Top of another loop point based on finding a linearly */
248  /* - independent species */
249  do {
250  /*
251  * Search the remaining part of the mole number vector, molNum
252  * for the largest remaining species. Return its identity.
253  * kk is the raw number. k is the orderVectorSpecies index.
254  */
255  kk = amax(DATA_PTR(molNum), 0, nspecies);
256  for (j = 0; j < nspecies; j++) {
257  if (orderVectorSpecies[j] == kk) {
258  k = j;
259  break;
260  }
261  }
262  if (j == nspecies) {
263  throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
264  }
265 
266  if (molNum[kk] == 0.0) {
267  *usedZeroedSpecies = true;
268  }
269  /*
270  * If the largest molNum is negative, then we are done.
271  */
272  if (molNum[kk] == USEDBEFORE) {
273  nComponents = jr;
274  nNonComponents = nspecies - nComponents;
275  goto L_END_LOOP;
276  }
277  /*
278  * Assign a small negative number to the component that we have
279  * just found, in order to take it out of further consideration.
280  */
281 #ifdef DEBUG_MODE
282  molSave = molNum[kk];
283 #endif
284  molNum[kk] = USEDBEFORE;
285 
286  /* *********************************************************** */
287  /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
288  /* *********************************************************** */
289  /*
290  * Modified Gram-Schmidt Method, p. 202 Dalquist
291  * QR factorization of a matrix without row pivoting.
292  */
293  jl = jr;
294  for (j = 0; j < ne; ++j) {
295  jj = orderVectorElements[j];
296  sm[j + jr*ne] = mphase->nAtoms(kk,jj);
297  }
298  if (jl > 0) {
299  /*
300  * Compute the coefficients of JA column of the
301  * the upper triangular R matrix, SS(J) = R_J_JR
302  * (this is slightly different than Dalquist)
303  * R_JA_JA = 1
304  */
305  for (j = 0; j < jl; ++j) {
306  ss[j] = 0.0;
307  for (i = 0; i < ne; ++i) {
308  ss[j] += sm[i + jr*ne] * sm[i + j*ne];
309  }
310  ss[j] /= sa[j];
311  }
312  /*
313  * Now make the new column, (*,JR), orthogonal to the
314  * previous columns
315  */
316  for (j = 0; j < jl; ++j) {
317  for (l = 0; l < ne; ++l) {
318  sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
319  }
320  }
321  }
322  /*
323  * Find the new length of the new column in Q.
324  * It will be used in the denominator in future row calcs.
325  */
326  sa[jr] = 0.0;
327  for (ml = 0; ml < ne; ++ml) {
328  tmp = sm[ml + jr*ne];
329  sa[jr] += tmp * tmp;
330  }
331  /* **************************************************** */
332  /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
333  /* **************************************************** */
334  if (sa[jr] < 1.0e-6) {
335  lindep = true;
336  } else {
337  lindep = false;
338  }
339  } while (lindep);
340  /* ****************************************** */
341  /* **** REARRANGE THE DATA ****************** */
342  /* ****************************************** */
343  if (jr != k) {
344 #ifdef DEBUG_MODE
345  if (BasisOptimize_print_lvl >= 1) {
346  kk = orderVectorSpecies[k];
347  sname = mphase->speciesName(kk);
348  writelogf(" --- %-12.12s", sname.c_str());
349  jj = orderVectorSpecies[jr];
350  ename = mphase->speciesName(jj);
351  writelogf("(%9.2g) replaces %-12.12s", molSave, ename.c_str());
352  writelogf("(%9.2g) as component %3d\n", molNum[jj], jr);
353  }
354 #endif
355  std::swap(orderVectorSpecies[jr], orderVectorSpecies[k]);
356  }
357  /* - entry point from up above */
358 L_END_LOOP:
359  ;
360  /*
361  * If we haven't found enough components, go back
362  * and find some more. (nc -1 is used below, because
363  * jr is counted from 0, via the C convention.
364  */
365  } while (jr < (nComponents-1));
366 
367 
368  if (! doFormRxn) {
369  return nComponents;
370  }
371 
372  /* ****************************************************** */
373  /* **** EVALUATE THE STOICHIOMETRY ********************** */
374  /* ****************************************************** */
375  /*
376  * Formulate the matrix problem for the stoichiometric
377  * coefficients. CX + B = 0
378  * C will be an nc x nc matrix made up of the formula
379  * vectors for the components. Each component's formula
380  * vector is a column. The rows are the elements.
381  * n rhs's will be solved for. Thus, B is an nc x n
382  * matrix.
383  *
384  * BIG PROBLEM 1/21/99:
385  *
386  * This algorithm makes the assumption that the
387  * first nc rows of the formula matrix aren't rank deficient.
388  * However, this might not be the case. For example, assume
389  * that the first element in FormulaMatrix[] is argon. Assume that
390  * no species in the matrix problem actually includes argon.
391  * Then, the first row in sm[], below will be indentically
392  * zero. bleh.
393  * What needs to be done is to perform a rearrangement
394  * of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
395  * that the first nc elements form in combination with the
396  * nc components create an invertible sm[]. not a small
397  * project, but very doable.
398  * An alternative would be to turn the matrix problem
399  * below into an ne x nc problem, and do QR elimination instead
400  * of Gauss-Jordon elimination.
401  * Note the rearrangement of elements need only be done once
402  * in the problem. It's actually very similar to the top of
403  * this program with ne being the species and nc being the
404  * elements!!
405  */
406  for (k = 0; k < nComponents; ++k) {
407  kk = orderVectorSpecies[k];
408  for (j = 0; j < nComponents; ++j) {
409  jj = orderVectorElements[j];
410  sm[j + k*ne] = mphase->nAtoms(kk, jj);
411  }
412  }
413 
414  for (i = 0; i < nNonComponents; ++i) {
415  k = nComponents + i;
416  kk = orderVectorSpecies[k];
417  for (j = 0; j < nComponents; ++j) {
418  jj = orderVectorElements[j];
419  formRxnMatrix[j + i * ne] = mphase->nAtoms(kk, jj);
420  }
421  }
422  /*
423  * Use Gauss-Jordon block elimination to calculate
424  * the reaction matrix
425  */
426  int ierr = mlequ(DATA_PTR(sm), ne, nComponents, DATA_PTR(formRxnMatrix), nNonComponents);
427  if (ierr == 1) {
428  writelog("ERROR: mlequ returned an error condition\n");
429  throw CanteraError("basopt", "mlequ returned an error condition");
430  }
431 
432 #ifdef DEBUG_MODE
433  if (Cantera::BasisOptimize_print_lvl >= 1) {
434  writelog(" ---\n");
435  writelogf(" --- Number of Components = %d\n", nComponents);
436  writelog(" --- Formula Matrix:\n");
437  writelog(" --- Components: ");
438  for (k = 0; k < nComponents; k++) {
439  kk = orderVectorSpecies[k];
440  writelogf(" %3d (%3d) ", k, kk);
441  }
442  writelog("\n --- Components Moles: ");
443  for (k = 0; k < nComponents; k++) {
444  kk = orderVectorSpecies[k];
445  writelogf("%-11.3g", molNumBase[kk]);
446  }
447  writelog("\n --- NonComponent | Moles | ");
448  for (i = 0; i < nComponents; i++) {
449  kk = orderVectorSpecies[i];
450  sname = mphase->speciesName(kk);
451  writelogf("%-11.10s", sname.c_str());
452  }
453  writelog("\n");
454 
455  for (i = 0; i < nNonComponents; i++) {
456  k = i + nComponents;
457  kk = orderVectorSpecies[k];
458  writelogf(" --- %3d (%3d) ", k, kk);
459  sname = mphase->speciesName(kk);
460  writelogf("%-10.10s", sname.c_str());
461  writelogf("|%10.3g|", molNumBase[kk]);
462  /*
463  * Print the negative of formRxnMatrix[]; it's easier to interpret.
464  */
465  for (j = 0; j < nComponents; j++) {
466  writelogf(" %6.2f", - formRxnMatrix[j + i * ne]);
467  }
468  writelog("\n");
469  }
470  writelog(" ");
471  for (i=0; i<77; i++) {
472  writelog("-");
473  }
474  writelog("\n");
475  }
476 #endif
477 
478  return nComponents;
479 } /* basopt() ************************************************************/
480 
481 
482 
483 #ifdef DEBUG_MODE
484 static void print_stringTrunc(const char* str, int space, int alignment)
485 
486 /***********************************************************************
487  * vcs_print_stringTrunc():
488  *
489  * Print a string within a given space limit. This routine
490  * limits the amount of the string that will be printed to a
491  * maximum of "space" characters.
492  *
493  * str = String -> must be null terminated.
494  * space = space limit for the printing.
495  * alignment = 0 centered
496  * 1 right aligned
497  * 2 left aligned
498  ***********************************************************************/
499 {
500  int i, ls=0, rs=0;
501  int len = strlen(str);
502  if ((len) >= space) {
503  for (i = 0; i < space; i++) {
504  writelogf("%c", str[i]);
505  }
506  } else {
507  if (alignment == 1) {
508  ls = space - len;
509  } else if (alignment == 2) {
510  rs = space - len;
511  } else {
512  ls = (space - len) / 2;
513  rs = space - len - ls;
514  }
515  if (ls != 0) {
516  for (i = 0; i < ls; i++) {
517  writelog(" ");
518  }
519  }
520  writelogf("%s", str);
521  if (rs != 0) {
522  for (i = 0; i < rs; i++) {
523  writelog(" ");
524  }
525  }
526  }
527 }
528 #endif
529 
530 /*
531  * Finds the location of the maximum component in a double vector
532  * INPUT
533  * x(*) - Vector to search
534  * j <= i < n : i is the range of indices to search in X(*)
535  *
536  * RETURN
537  * return index of the greatest value on X(*) searched
538  */
539 static size_t amax(double* x, size_t j, size_t n)
540 {
541  size_t largest = j;
542  double big = x[j];
543  for (size_t i = j + 1; i < n; ++i) {
544  if (x[i] > big) {
545  largest = i;
546  big = x[i];
547  }
548  }
549  return largest;
550 }
551 
552 /*
553  * vcs_mlequ:
554  *
555  * Invert an nxn matrix and solve m rhs's
556  *
557  * Solve C X + B = 0;
558  *
559  * This routine uses Gauss elimination and is optimized for the solution
560  * of lots of rhs's.
561  * A crude form of row pivoting is used here.
562  *
563  *
564  * c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number
565  * j = column number
566  * b[i+j*idem] = b_i_j = vectors of rhs's: i = row number
567  * j = column number
568  * (each column is a new rhs)
569  * n = number of rows and columns in the matrix
570  * m = number of rhs to be solved for
571  * idem = first dimension in the calling routine
572  * idem >= n must be true
573  *
574  * Return Value
575  * 1 : Matrix is singular
576  * 0 : solution is OK
577  *
578  * The solution is returned in the matrix b.
579  */
580 static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m)
581 {
582  size_t i, j, k, l;
583  double R;
584 
585  /*
586  * Loop over the rows
587  * -> At the end of each loop, the only nonzero entry in the column
588  * will be on the diagonal. We can therfore just invert the
589  * diagonal at the end of the program to solve the equation system.
590  */
591  for (i = 0; i < n; ++i) {
592  if (c[i + i * idem] == 0.0) {
593  /*
594  * Do a simple form of row pivoting to find a non-zero pivot
595  */
596  for (k = i + 1; k < n; ++k) {
597  if (c[k + i * idem] != 0.0) {
598  goto FOUND_PIVOT;
599  }
600  }
601 #ifdef DEBUG_MODE
602  writelogf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i);
603 #endif
604  return 1;
605 FOUND_PIVOT:
606  ;
607  for (j = 0; j < n; ++j) {
608  c[i + j * idem] += c[k + j * idem];
609  }
610  for (j = 0; j < m; ++j) {
611  b[i + j * idem] += b[k + j * idem];
612  }
613  }
614 
615  for (l = 0; l < n; ++l) {
616  if (l != i && c[l + i * idem] != 0.0) {
617  R = c[l + i * idem] / c[i + i * idem];
618  c[l + i * idem] = 0.0;
619  for (j = i+1; j < n; ++j) {
620  c[l + j * idem] -= c[i + j * idem] * R;
621  }
622  for (j = 0; j < m; ++j) {
623  b[l + j * idem] -= b[i + j * idem] * R;
624  }
625  }
626  }
627  }
628  /*
629  * The negative in the last expression is due to the form of B upon
630  * input
631  */
632  for (i = 0; i < n; ++i) {
633  for (j = 0; j < m; ++j) {
634  b[i + j * idem] = -b[i + j * idem] / c[i + i*idem];
635  }
636  }
637  return 0;
638 } /* mlequ() *************************************************************/
639 
640 
641 /*
642  *
643  * ElemRearrange:
644  *
645  * This subroutine handles the rearrangement of the constraint
646  * equations represented by the Formula Matrix. Rearrangement is only
647  * necessary when the number of components is less than the number of
648  * elements. For this case, some constraints can never be satisfied
649  * exactly, because the range space represented by the Formula
650  * Matrix of the components can't span the extra space. These
651  * constraints, which are out of the range space of the component
652  * Formula matrix entries, are migrated to the back of the Formula
653  * matrix.
654  *
655  * A prototypical example is an extra element column in
656  * FormulaMatrix[],
657  * which is identically zero. For example, let's say that argon is
658  * has an element column in FormulaMatrix[], but no species in the
659  * mechanism
660  * actually contains argon. Then, nc < ne. Unless the entry for
661  * desired element abundance vector for Ar is zero, then this
662  * element abundance constraint can never be satisfied. The
663  * constraint vector is not in the range space of the formula
664  * matrix.
665  * Also, without perturbation
666  * of FormulaMatrix[], BasisOptimize[] would produce a zero pivot
667  * because the matrix
668  * would be singular (unless the argon element column was already the
669  * last column of FormulaMatrix[].
670  * This routine borrows heavily from BasisOptimize algorithm. It
671  * finds nc constraints which span the range space of the Component
672  * Formula matrix, and assigns them as the first nc components in the
673  * formula matrix. This guarantees that BasisOptimize has a
674  * nonsingular matrix to invert.
675  */
676 size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
677  MultiPhase* mphase,
678  std::vector<size_t>& orderVectorSpecies,
679  std::vector<size_t>& orderVectorElements)
680 {
681 
682  size_t j, k, l, i, jl, ml, jr, ielem, jj, kk=0;
683 
684  bool lindep = false;
685  size_t nelements = mphase->nElements();
686  std::string ename;
687  /*
688  * Get the total number of species in the multiphase object
689  */
690  size_t nspecies = mphase->nSpecies();
691 
692  double test = -1.0E10;
693 #ifdef DEBUG_MODE
694  if (BasisOptimize_print_lvl > 0) {
695  writelog(" ");
696  for (i=0; i<77; i++) {
697  writelog("-");
698  }
699  writelog("\n");
700  writelog(" --- Subroutine ElemRearrange() called to ");
701  writelog("check stoich. coefficient matrix\n");
702  writelog(" --- and to rearrange the element ordering once\n");
703  }
704 #endif
705 
706  /*
707  * Perhaps, initialize the element ordering
708  */
709  if (orderVectorElements.size() < nelements) {
710  orderVectorElements.resize(nelements);
711  for (j = 0; j < nelements; j++) {
712  orderVectorElements[j] = j;
713  }
714  }
715 
716  /*
717  * Perhaps, initialize the species ordering. However, this is
718  * dangerous, as this ordering is assumed to yield the
719  * component species for the problem
720  */
721  if (orderVectorSpecies.size() != nspecies) {
722  orderVectorSpecies.resize(nspecies);
723  for (k = 0; k < nspecies; k++) {
724  orderVectorSpecies[k] = k;
725  }
726  }
727 
728  /*
729  * If the elementAbundances aren't input, just create a fake one
730  * based on summing the column of the stoich matrix.
731  * This will force elements with zero species to the
732  * end of the element ordering.
733  */
734  vector_fp eAbund(nelements,0.0);
735  if (elementAbundances.size() != nelements) {
736  for (j = 0; j < nelements; j++) {
737  eAbund[j] = 0.0;
738  for (k = 0; k < nspecies; k++) {
739  eAbund[j] += fabs(mphase->nAtoms(k, j));
740  }
741  }
742  } else {
743  copy(elementAbundances.begin(), elementAbundances.end(),
744  eAbund.begin());
745  }
746 
747  vector_fp sa(nelements,0.0);
748  vector_fp ss(nelements,0.0);
749  vector_fp sm(nelements*nelements,0.0);
750 
751  /*
752  * Top of a loop of some sort based on the index JR. JR is the
753  * current number independent elements found.
754  */
755  jr = npos;
756  do {
757  ++jr;
758  /*
759  * Top of another loop point based on finding a linearly
760  * independent element
761  */
762  do {
763  /*
764  * Search the element vector. We first locate elements that
765  * are present in any amount. Then, we locate elements that
766  * are not present in any amount.
767  * Return its identity in K.
768  */
769  k = nelements;
770  for (ielem = jr; ielem < nelements; ielem++) {
771  kk = orderVectorElements[ielem];
772  if (eAbund[kk] != test && eAbund[kk] > 0.0) {
773  k = ielem;
774  break;
775  }
776  }
777  for (ielem = jr; ielem < nelements; ielem++) {
778  kk = orderVectorElements[ielem];
779  if (eAbund[kk] != test) {
780  k = ielem;
781  break;
782  }
783  }
784 
785  if (k == nelements) {
786  // When we are here, there is an error usually.
787  // We haven't found the number of elements necessary.
788  // This is signalled by returning jr != nComponents.
789 #ifdef DEBUG_MODE
790  if (BasisOptimize_print_lvl > 0) {
791  writelogf("Error exit: returning with nComponents = %d\n", jr);
792  }
793 #endif
794  return jr;
795  }
796 
797  /*
798  * Assign a large negative number to the element that we have
799  * just found, in order to take it out of further consideration.
800  */
801  eAbund[kk] = test;
802 
803  /* *********************************************************** */
804  /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */
805  /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */
806  /* *********************************************************** */
807  /*
808  * Modified Gram-Schmidt Method, p. 202 Dalquist
809  * QR factorization of a matrix without row pivoting.
810  */
811  jl = jr;
812  /*
813  * Fill in the row for the current element, k, under consideration
814  * The row will contain the Formula matrix value for that element
815  * with respect to the vector of component species.
816  * (note j and k indices are flipped compared to the previous routine)
817  */
818  for (j = 0; j < nComponents; ++j) {
819  jj = orderVectorSpecies[j];
820  kk = orderVectorElements[k];
821  sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
822  }
823  if (jl > 0) {
824  /*
825  * Compute the coefficients of JA column of the
826  * the upper triangular R matrix, SS(J) = R_J_JR
827  * (this is slightly different than Dalquist)
828  * R_JA_JA = 1
829  */
830  for (j = 0; j < jl; ++j) {
831  ss[j] = 0.0;
832  for (i = 0; i < nComponents; ++i) {
833  ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
834  }
835  ss[j] /= sa[j];
836  }
837  /*
838  * Now make the new column, (*,JR), orthogonal to the
839  * previous columns
840  */
841  for (j = 0; j < jl; ++j) {
842  for (l = 0; l < nComponents; ++l) {
843  sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
844  }
845  }
846  }
847 
848  /*
849  * Find the new length of the new column in Q.
850  * It will be used in the denominator in future row calcs.
851  */
852  sa[jr] = 0.0;
853  for (ml = 0; ml < nComponents; ++ml) {
854  double tmp = sm[ml + jr*nComponents];
855  sa[jr] += tmp * tmp;
856  }
857  /* **************************************************** */
858  /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */
859  /* **************************************************** */
860  if (sa[jr] < 1.0e-6) {
861  lindep = true;
862  } else {
863  lindep = false;
864  }
865  } while (lindep);
866  /* ****************************************** */
867  /* **** REARRANGE THE DATA ****************** */
868  /* ****************************************** */
869  if (jr != k) {
870 #ifdef DEBUG_MODE
871  if (BasisOptimize_print_lvl > 0) {
872  kk = orderVectorElements[k];
873  ename = mphase->elementName(kk);
874  writelog(" --- ");
875  writelogf("%-2.2s", ename.c_str());
876  writelog("replaces ");
877  kk = orderVectorElements[jr];
878  ename = mphase->elementName(kk);
879  writelogf("%-2.2s", ename.c_str());
880  writelogf(" as element %3d\n", jr);
881  }
882 #endif
883  std::swap(orderVectorElements[jr], orderVectorElements[k]);
884  }
885 
886  /*
887  * If we haven't found enough components, go back
888  * and find some more. (nc -1 is used below, because
889  * jr is counted from 0, via the C convention.
890  */
891  } while (jr < (nComponents-1));
892  return nComponents;
893 } /* vcs_elem_rearrange() ****************************************************/
894 
Generated by   doxygen 1.8.2