BasisOptimize.cpp

Go to the documentation of this file.
00001 /**
00002  * @file BasisOptimize.cpp
00003  *     Functions which calculation optimized basis of the 
00004  *     stoichiometric coefficient matrix (see /ref equil functions)
00005  */
00006 
00007 /*
00008  * $Id: BasisOptimize.cpp 368 2010-01-04 00:46:26Z hkmoffa $
00009  */
00010 
00011 #include "ct_defs.h"
00012 #include "global.h"
00013 #include "ThermoPhase.h"
00014 #include "MultiPhase.h"
00015 
00016 #include <string.h>
00017 
00018 using namespace Cantera;
00019 using namespace std;
00020 
00021 #ifdef DEBUG_MODE
00022 namespace Cantera {
00023   int BasisOptimize_print_lvl = 0;
00024 }
00025 static void print_stringTrunc(const char *str, int space, int alignment);
00026 #endif
00027 
00028 static int amax(double *x, int j, int n);
00029 static void switch_pos(vector_int &orderVector, int jr, int kspec);
00030 static int mlequ(double *c, int idem, int n, double *b, int m);
00031 
00032 //@{
00033 #ifndef MIN
00034 #define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
00035 #endif
00036 //@}
00037 
00038 /*
00039  * Choose the optimum basis for the calculations. This is done by 
00040  * choosing the species with the largest mole fraction 
00041  * not currently a linear combination of the previous components. 
00042  * Then, calculate the stoichiometric coefficient matrix for that 
00043  * basis. 
00044  *
00045  * Calculates the identity of the component species in the mechanism. 
00046  * Rearranges the solution data to put the component data at the 
00047  * front of the species list. 
00048  *
00049  * Then, calculates SC(J,I) the formation reactions for all noncomponent 
00050  * species in the mechanism. 
00051  *
00052  * Input 
00053  * --------- 
00054  * mphase          Pointer to the multiphase object. Contains the 
00055  *                 species mole fractions, which are used to pick the
00056  *                 current optimal species component basis.
00057  * orderVectorElement
00058  *                 Order vector for the elements. The element rows
00059  *                 in the formula matrix are
00060  *                 rearranged according to this vector.
00061  * orderVectorSpecies
00062  *                 Order vector for the species. The species are
00063  *                 rearranged according to this formula. The first
00064  *                 nCompoments of this vector contain the calculated
00065  *                 species components on exit.
00066  * doFormRxn       If true, the routine calculates the formation
00067  *                 reaction matrix based on the calculated 
00068  *                 component species. If false, this step is skipped.
00069  * 
00070  * Output 
00071  * --------- 
00072  * usedZeroedSpecies = If true, then a species with a zero concentration
00073  *                     was used as a component. The problem may be
00074  *                     converged.
00075  * formRxnMatrix 
00076  *
00077  * Return
00078  * --------------
00079  * returns the number of components.
00080  *
00081  *
00082  */
00083 int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn,
00084                MultiPhase *mphase, vector_int & orderVectorSpecies,
00085                vector_int & orderVectorElements, 
00086                vector_fp & formRxnMatrix) {
00087 
00088   int  j, jj, k=0, kk, l, i, jl, ml;
00089   bool lindep;
00090   std::string ename;
00091   std::string sname;
00092   /*
00093    * Get the total number of elements defined in the multiphase object
00094    */
00095   int ne = mphase->nElements();
00096   /*
00097    * Get the total number of species in the multiphase object
00098    */
00099   int nspecies = mphase->nSpecies();
00100   doublereal tmp;
00101   doublereal const USEDBEFORE = -1;
00102  
00103   /*
00104    * Perhaps, initialize the element ordering
00105    */
00106   if ((int) orderVectorElements.size() < ne) {
00107     orderVectorElements.resize(ne);
00108     for (j = 0; j < ne; j++) {
00109       orderVectorElements[j] = j;
00110     }
00111   }
00112 
00113   /*
00114    * Perhaps, initialize the species ordering
00115    */
00116   if ((int) orderVectorSpecies.size() != nspecies) {
00117     orderVectorSpecies.resize(nspecies);
00118     for (k = 0; k < nspecies; k++) {
00119       orderVectorSpecies[k] = k;
00120     }
00121   }
00122  
00123 #ifdef DEBUG_MODE
00124   double molSave = 0.0;
00125   if (BasisOptimize_print_lvl >= 1) {
00126     writelog("   "); for(i=0; i<77; i++) writelog("-"); writelog("\n");
00127     writelog("   --- Subroutine BASOPT called to ");
00128     writelog("calculate the number of components and ");
00129     writelog("evaluate the formation matrix\n");
00130     if (BasisOptimize_print_lvl > 0) {
00131       writelog("   ---\n");
00132       
00133       writelog("   ---      Formula Matrix used in BASOPT calculation\n");
00134       writelog("   ---      Species | Order | ");
00135       for (j = 0; j < ne; j++) {
00136     jj = orderVectorElements[j];
00137     writelog(" ");
00138     ename = mphase->elementName(jj);
00139     print_stringTrunc(ename.c_str(), 4, 1);
00140     writelogf("(%1d)", j); 
00141       }
00142       writelog("\n");
00143       for (k = 0; k < nspecies; k++) {
00144     kk = orderVectorSpecies[k];
00145     writelog("   --- ");
00146     sname = mphase->speciesName(kk);
00147     print_stringTrunc(sname.c_str(), 11, 1);
00148     writelogf(" |   %4d |", k); 
00149     for (j = 0; j < ne; j++) {
00150       jj = orderVectorElements[j]; 
00151       double num = mphase->nAtoms(kk,jj);
00152       writelogf("%6.1g  ", num); 
00153     }
00154     writelog("\n");
00155       }
00156       writelog("   --- \n");
00157     }
00158   }
00159 #endif
00160    
00161   /*
00162    *  Calculate the maximum value of the number of components possible
00163    *     It's equal to the minimum of the number of elements and the
00164    *     number of total species.
00165    */
00166   int nComponents = MIN(ne, nspecies);
00167   int nNonComponents = nspecies - nComponents;
00168   /*
00169    * Set this return variable to false
00170    */
00171   *usedZeroedSpecies = false;
00172 
00173   /*
00174    * Create an array of mole numbers
00175    */
00176   vector_fp molNum(nspecies,0.0);
00177   mphase->getMoles(DATA_PTR(molNum));
00178 
00179   /*
00180    * Other workspace
00181    */
00182   vector_fp sm(ne*ne, 0.0);
00183   vector_fp ss(ne, 0.0);
00184   vector_fp sa(ne, 0.0);
00185   if ((int) formRxnMatrix.size() < nspecies*ne) {
00186     formRxnMatrix.resize(nspecies*ne, 0.0);
00187   }
00188 
00189 #ifdef DEBUG_MODE
00190   /*
00191    * For debugging purposes keep an unmodified copy of the array.
00192    */
00193   vector_fp molNumBase(molNum);
00194 #endif
00195   
00196 
00197   int jr = -1;
00198   /*
00199    *   Top of a loop of some sort based on the index JR. JR is the 
00200    *   current number of component species found. 
00201    */
00202   do {
00203     ++jr;
00204     /* - Top of another loop point based on finding a linearly */
00205     /* - independent species */
00206     do {
00207       /*
00208        *    Search the remaining part of the mole number vector, molNum
00209        *    for the largest remaining species. Return its identity. 
00210        *    kk is the raw number. k is the orderVectorSpecies index.
00211        */
00212       kk = amax(DATA_PTR(molNum), 0, nspecies);
00213       for (j = 0; j < nspecies; j++) {
00214     if (orderVectorSpecies[j] == kk) {
00215       k = j;
00216       break;
00217     }
00218       }
00219       if (j == nspecies) {
00220     throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
00221       }
00222 
00223       if (molNum[kk] == 0.0) *usedZeroedSpecies = true;
00224       /*
00225        * If the largest molNum is negative, then we are done.
00226        */
00227       if (molNum[kk] == USEDBEFORE) {
00228     nComponents = jr;
00229     nNonComponents = nspecies - nComponents;
00230     goto L_END_LOOP;
00231       }
00232       /*
00233        *  Assign a small negative number to the component that we have
00234        *  just found, in order to take it out of further consideration.
00235        */
00236 #ifdef DEBUG_MODE
00237       molSave = molNum[kk];
00238 #endif
00239       molNum[kk] = USEDBEFORE;
00240 
00241       /* *********************************************************** */
00242       /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
00243       /* *********************************************************** */
00244       /*    
00245        *          Modified Gram-Schmidt Method, p. 202 Dalquist 
00246        *          QR factorization of a matrix without row pivoting. 
00247        */
00248       jl = jr;
00249       for (j = 0; j < ne; ++j) {
00250     jj = orderVectorElements[j];
00251     sm[j + jr*ne] = mphase->nAtoms(kk,jj);
00252       }
00253       if (jl > 0) {
00254     /*
00255      *         Compute the coefficients of JA column of the 
00256      *         the upper triangular R matrix, SS(J) = R_J_JR 
00257      *         (this is slightly different than Dalquist) 
00258      *         R_JA_JA = 1 
00259      */
00260     for (j = 0; j < jl; ++j) {
00261       ss[j] = 0.0;
00262       for (i = 0; i < ne; ++i) {
00263         ss[j] += sm[i + jr*ne] * sm[i + j*ne];
00264       }
00265       ss[j] /= sa[j];
00266     }
00267     /* 
00268      *     Now make the new column, (*,JR), orthogonal to the 
00269      *     previous columns
00270      */
00271     for (j = 0; j < jl; ++j) {
00272       for (l = 0; l < ne; ++l) {
00273         sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
00274       }
00275     }
00276       }
00277       /*
00278        *        Find the new length of the new column in Q. 
00279        *        It will be used in the denominator in future row calcs. 
00280        */
00281       sa[jr] = 0.0;
00282       for (ml = 0; ml < ne; ++ml) {
00283     tmp = sm[ml + jr*ne];
00284     sa[jr] += tmp * tmp;
00285       }
00286       /* **************************************************** */
00287       /* **** IF NORM OF NEW ROW  .LT. 1E-3 REJECT ********** */
00288       /* **************************************************** */
00289       if (sa[jr] < 1.0e-6)  lindep = true;
00290       else                  lindep = false;
00291     } while(lindep);
00292     /* ****************************************** */
00293     /* **** REARRANGE THE DATA ****************** */
00294     /* ****************************************** */
00295     if (jr != k) {
00296 #ifdef DEBUG_MODE
00297       if (BasisOptimize_print_lvl >= 1) {
00298     kk = orderVectorSpecies[k];
00299     sname = mphase->speciesName(kk);
00300     writelogf("   ---   %-12.12s", sname.c_str()); 
00301     jj = orderVectorSpecies[jr];
00302     ename = mphase->speciesName(jj);
00303     writelogf("(%9.2g) replaces %-12.12s", molSave, ename.c_str());
00304     writelogf("(%9.2g) as component %3d\n", molNum[jj], jr);
00305       }
00306 #endif
00307       switch_pos(orderVectorSpecies, jr, k);
00308     }
00309     /* - entry point from up above */
00310   L_END_LOOP: ;
00311     /*
00312      *      If we haven't found enough components, go back 
00313      *      and find some more. (nc -1 is used below, because
00314      *      jr is counted from 0, via the C convention.
00315      */
00316   } while (jr < (nComponents-1));
00317    
00318 
00319  if (! doFormRxn) return nComponents;
00320 
00321   /* ****************************************************** */
00322   /* **** EVALUATE THE STOICHIOMETRY ********************** */
00323   /* ****************************************************** */
00324   /*
00325    *  Formulate the matrix problem for the stoichiometric
00326    *  coefficients. CX + B = 0
00327    *      C will be an nc x nc matrix made up of the formula 
00328    * vectors for the components. Each component's formular
00329    * vector is a column. The rows are the elements.
00330    *      n rhs's will be solved for. Thus, B is an nc x n
00331    * matrix. 
00332    *
00333    * BIG PROBLEM 1/21/99:
00334    *
00335    *    This algorithm makes the assumption that the
00336    * first nc rows of the formula matrix aren't rank deficient.
00337    * However, this might not be the case. For example, assume
00338    * that the first element in FormulaMatrix[] is argon. Assume that
00339    * no species in the matrix problem actually includes argon.
00340    * Then, the first row in sm[], below will be indentically
00341    * zero. bleh. 
00342    *    What needs to be done is to perform a rearrangement
00343    * of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
00344    * that the first nc elements form in combination with the
00345    * nc components create an invertible sm[]. not a small
00346    * project, but very doable.
00347    *    An alternative would be to turn the matrix problem
00348    * below into an ne x nc problem, and do QR elimination instead
00349    * of Gauss-Jordon elimination.
00350    *    Note the rearrangement of elements need only be done once
00351    * in the problem. It's actually very similar to the top of 
00352    * this program with ne being the species and nc being the
00353    * elements!!
00354    */
00355   for (k = 0; k < nComponents; ++k) {
00356     kk = orderVectorSpecies[k];
00357     for (j = 0; j < nComponents; ++j) {
00358       jj = orderVectorElements[j];
00359       sm[j + k*ne] = mphase->nAtoms(kk, jj);
00360     }
00361   }
00362 
00363   for (i = 0; i < nNonComponents; ++i) {
00364     k = nComponents + i;
00365     kk = orderVectorSpecies[k];
00366     for (j = 0; j < nComponents; ++j) {
00367       jj = orderVectorElements[j];
00368       formRxnMatrix[j + i * ne] = mphase->nAtoms(kk, jj);
00369     }
00370   }
00371   /*
00372    *     Use Gauss-Jordon block elimination to calculate
00373    *     the reaction matrix 
00374    */
00375   j = mlequ(DATA_PTR(sm), ne, nComponents, DATA_PTR(formRxnMatrix), nNonComponents);
00376   if (j == 1) {
00377     writelog("ERROR: mlequ returned an error condition\n");
00378     throw CanteraError("basopt", "mlequ returned an error condition");
00379   }
00380     
00381 #ifdef DEBUG_MODE
00382   if (Cantera::BasisOptimize_print_lvl >= 1) {
00383     writelog("   ---\n");
00384     writelogf("   ---  Number of Components = %d\n", nComponents);
00385     writelog("   ---  Formula Matrix:\n");
00386     writelog("   ---                      Components:    ");
00387     for (k = 0; k < nComponents; k++) {
00388       kk = orderVectorSpecies[k];
00389       writelogf(" %3d (%3d) ", k, kk);
00390     }
00391     writelog("\n   ---                Components Moles:       ");
00392     for (k = 0; k < nComponents; k++) {
00393       kk = orderVectorSpecies[k];
00394       writelogf("%-11.3g", molNumBase[kk]); 
00395     }
00396     writelog("\n   ---        NonComponent |   Moles  |       ");
00397     for (i = 0; i < nComponents; i++) {
00398       kk = orderVectorSpecies[i];
00399       sname = mphase->speciesName(kk);
00400       writelogf("%-11.10s", sname.c_str());
00401     }
00402     writelog("\n");
00403   
00404     for (i = 0; i < nNonComponents; i++) {
00405       k = i + nComponents;
00406       kk = orderVectorSpecies[k];
00407       writelogf("   --- %3d (%3d) ", k, kk);
00408       sname = mphase->speciesName(kk);
00409       writelogf("%-10.10s", sname.c_str()); 
00410       writelogf("|%10.3g|", molNumBase[kk]); 
00411       /*
00412        * Print the negative of formRxnMatrix[]; it's easier to interpret.
00413        */
00414       for (j = 0; j < nComponents; j++) {
00415     writelogf("     %6.2f", - formRxnMatrix[j + i * ne]);
00416       }
00417       writelog("\n");
00418     }
00419     writelog("   "); for (i=0; i<77; i++) writelog("-"); writelog("\n");
00420   }
00421 #endif
00422 
00423   return nComponents;
00424 } /* basopt() ************************************************************/
00425 
00426 
00427 
00428 #ifdef DEBUG_MODE
00429 static void print_stringTrunc(const char *str, int space, int alignment)
00430 
00431    /***********************************************************************
00432     *  vcs_print_stringTrunc():
00433     *
00434     *     Print a string within a given space limit. This routine
00435     *     limits the amount of the string that will be printed to a
00436     *     maximum of "space" characters.
00437     *
00438     *     str = String -> must be null terminated.
00439     *     space = space limit for the printing.
00440     *     alignment = 0 centered
00441     *           1 right aligned
00442     *           2 left aligned
00443     ***********************************************************************/
00444 {
00445   int i, ls=0, rs=0;
00446   int len = strlen(str);
00447   if ((len) >= space) {
00448     for (i = 0; i < space; i++) {
00449       writelogf("%c", str[i]);
00450     }
00451   } else {
00452     if (alignment == 1) {
00453       ls = space - len;
00454     } else if (alignment == 2) {
00455       rs = space - len;
00456     } else {
00457       ls = (space - len) / 2;
00458       rs = space - len - ls;
00459     }
00460     if (ls != 0) {
00461       for (i = 0; i < ls; i++) writelog(" ");
00462     }
00463     writelogf("%s", str);
00464     if (rs != 0) {
00465       for (i = 0; i < rs; i++) writelog(" ");
00466     }
00467   }
00468 }
00469 #endif
00470 
00471 /*
00472  * Finds the location of the maximum component in a double vector
00473  * INPUT
00474  *    x(*) - Vector to search
00475  *    j <= i < n     : i is the range of indecises to search in X(*)
00476  *
00477  * RETURN
00478  *    return index of the greatest value on X(*) searched
00479  */
00480 static int amax(double *x, int j, int n) {
00481   int i;
00482   int largest = j;
00483   double big = x[j];
00484   for (i = j + 1; i < n; ++i) {
00485     if (x[i] > big) {
00486       largest = i;
00487       big = x[i];
00488     }
00489   }
00490   return largest;
00491 }
00492 
00493 
00494  static void switch_pos(vector_int &orderVector, int jr, int kspec) {
00495    int kcurr = orderVector[jr];
00496    orderVector[jr] = orderVector[kspec];
00497    orderVector[kspec] = kcurr;
00498  }
00499 
00500  /*
00501   * vcs_mlequ:
00502   *
00503   *  Invert an nxn matrix and solve m rhs's
00504   *
00505   *    Solve         C X + B = 0;
00506   *
00507   * This routine uses Gauss elimination and is optimized for the solution
00508   * of lots of rhs's.
00509   * A crude form of row pivoting is used here.
00510   *
00511   *
00512   * c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number
00513   *                                              j = column number
00514   * b[i+j*idem] = b_i_j = vectors of rhs's:      i = row number
00515   *                                              j = column number
00516   *            (each column is a new rhs)
00517   * n = number of rows and columns in the matrix
00518   * m = number of rhs to be solved for
00519   * idem = first dimension in the calling routine
00520   *        idem >= n must be true
00521   *
00522   * Return Value
00523   *      1 : Matrix is singluar
00524   *      0 : solution is OK
00525   *
00526   *      The solution is returned in the matrix b.
00527   */
00528  static int mlequ(double *c, int idem, int n, double *b, int m) {
00529    int i, j, k, l;
00530    double R;
00531 
00532    /*
00533     * Loop over the rows
00534     *    -> At the end of each loop, the only nonzero entry in the column
00535     *       will be on the diagonal. We can therfore just invert the
00536     *       diagonal at the end of the program to solve the equation system.
00537     */
00538    for (i = 0; i < n; ++i) {
00539      if (c[i + i * idem] == 0.0) {
00540        /*
00541     *   Do a simple form of row pivoting to find a non-zero pivot
00542     */
00543        for (k = i + 1; k < n; ++k) {
00544      if (c[k + i * idem] != 0.0) goto FOUND_PIVOT;
00545        }
00546 #ifdef DEBUG_MODE
00547        writelogf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i); 
00548 #endif
00549        return 1;
00550      FOUND_PIVOT: ;
00551        for (j = 0; j < n; ++j) c[i + j * idem] += c[k + j * idem];
00552        for (j = 0; j < m; ++j) b[i + j * idem] += b[k + j * idem];
00553      }
00554 
00555      for (l = 0; l < n; ++l) {
00556        if (l != i && c[l + i * idem] != 0.0) {
00557      R = c[l + i * idem] / c[i + i * idem];
00558      c[l + i * idem] = 0.0;
00559      for (j = i+1; j < n; ++j) c[l + j * idem] -= c[i + j * idem] * R;
00560      for (j = 0; j < m; ++j)   b[l + j * idem] -= b[i + j * idem] * R;
00561        }
00562      }
00563    }
00564    /*
00565     *  The negative in the last expression is due to the form of B upon
00566     *  input
00567     */
00568    for (i = 0; i < n; ++i) {
00569      for (j = 0; j < m; ++j)
00570        b[i + j * idem] = -b[i + j * idem] / c[i + i*idem];
00571    }
00572    return 0;
00573  } /* mlequ() *************************************************************/
00574 
00575 
00576 /*
00577  *
00578  * ElemRearrange:
00579  *
00580  *    This subroutine handles the rearrangement of the constraint
00581  *    equations represented by the Formula Matrix. Rearrangement is only
00582  *    necessary when the number of components is less than the number of
00583  *    elements. For this case, some constraints can never be satisfied 
00584  *    exactly, because the range space represented by the Formula
00585  *    Matrix of the components can't span the extra space. These 
00586  *    constraints, which are out of the range space of the component
00587  *    Formula matrix entries, are migrated to the back of the Formula
00588  *    matrix.
00589  *
00590  *    A prototypical example is an extra element column in 
00591  *    FormulaMatrix[], 
00592  *    which is identically zero. For example, let's say that argon is
00593  *    has an element column in FormulaMatrix[], but no species in the 
00594  *    mechanism
00595  *    actually contains argon. Then, nc < ne. Unless the entry for
00596  *    desired elementabundance vector for Ar is zero, then this
00597  *    element abundance constraint can never be satisfied. The 
00598  *    constraint vector is not in the range space of the formula
00599  *    matrix.
00600  *    Also, without perturbation
00601  *    of FormulaMatrix[], BasisOptimize[] would produce a zero pivot 
00602  *    because the matrix
00603  *    would be singular (unless the argon element column was already the
00604  *    last column of  FormulaMatrix[]. 
00605  *       This routine borrows heavily from BasisOptimize algorithm. It 
00606  *    finds nc constraints which span the range space of the Component
00607  *    Formula matrix, and assigns them as the first nc components in the
00608  *    formular matrix. This guarrantees that BasisOptimize has a
00609  *    nonsingular matrix to invert.
00610  */
00611 int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances,
00612                MultiPhase *mphase, 
00613                vector_int & orderVectorSpecies,
00614                vector_int & orderVectorElements) {
00615  
00616   int  j, k, l, i, jl, ml, jr, ielem, jj, kk=0;
00617  
00618   bool lindep = false;
00619   int nelements = mphase->nElements();
00620   std::string ename;
00621   /*
00622    * Get the total number of species in the multiphase object
00623    */
00624   int nspecies = mphase->nSpecies();
00625 
00626   double test = -1.0E10;
00627 #ifdef DEBUG_MODE
00628   if (BasisOptimize_print_lvl > 0) {
00629     writelog("   "); for(i=0; i<77; i++) writelog("-"); writelog("\n");
00630     writelog("   --- Subroutine ElemRearrange() called to ");
00631     writelog("check stoich. coefficent matrix\n");
00632     writelog("   ---    and to rearrange the element ordering once\n");
00633   }
00634 #endif
00635 
00636   /*
00637    * Perhaps, initialize the element ordering
00638    */
00639   if ((int) orderVectorElements.size() < nelements) {
00640     orderVectorElements.resize(nelements);
00641     for (j = 0; j < nelements; j++) {
00642       orderVectorElements[j] = j;
00643     }
00644   }
00645 
00646   /*
00647    * Perhaps, initialize the species ordering. However, this is 
00648    * dangerous, as this ordering is assumed to yield the
00649    * component species for the problem
00650    */
00651   if ((int) orderVectorSpecies.size() != nspecies) {
00652     orderVectorSpecies.resize(nspecies);
00653     for (k = 0; k < nspecies; k++) {
00654       orderVectorSpecies[k] = k;
00655     }
00656   }
00657 
00658   /*
00659    * If the elementAbundances aren't input, just create a fake one
00660    * based on summing the column of the stoich matrix.
00661    * This will force elements with zero species to the
00662    * end of the element ordering.
00663    */
00664   vector_fp eAbund(nelements,0.0);
00665   if ((int) elementAbundances.size() != nelements) {
00666     for (j = 0; j < nelements; j++) {
00667       eAbund[j] = 0.0;
00668       for (k = 0; k < nspecies; k++) {
00669     eAbund[j] += fabs(mphase->nAtoms(k, j));
00670       }
00671     }
00672   } else {
00673     copy(elementAbundances.begin(), elementAbundances.end(), 
00674      eAbund.begin());
00675   }
00676 
00677   vector_fp sa(nelements,0.0);
00678   vector_fp ss(nelements,0.0);
00679   vector_fp sm(nelements*nelements,0.0);
00680    
00681   /*
00682    *        Top of a loop of some sort based on the index JR. JR is the 
00683    *       current number independent elements found. 
00684    */
00685   jr = -1;
00686   do {
00687     ++jr;
00688     /* 
00689      *     Top of another loop point based on finding a linearly 
00690      *     independent element
00691      */
00692     do {
00693       /*
00694        *    Search the element vector. We first locate elements that
00695        *    are present in any amount. Then, we locate elements that
00696        *    are not present in any amount.
00697        *    Return its identity in K. 
00698        */
00699       k = nelements;
00700       for (ielem = jr; ielem < nelements; ielem++) {
00701     kk = orderVectorElements[ielem];
00702     if (eAbund[kk] != test && eAbund[kk] > 0.0) {
00703       k = ielem;
00704       break;
00705     }
00706       }
00707       for (ielem = jr; ielem < nelements; ielem++) {
00708     kk = orderVectorElements[ielem];
00709     if (eAbund[kk] != test) {
00710       k = ielem;
00711       break;
00712     }
00713       }
00714 
00715       if (k == nelements) {
00716     // When we are here, there is an error usually.
00717     // We haven't found the number of elements necessary.
00718     // This is signalled by returning jr != nComponents.
00719 #ifdef DEBUG_MODE
00720       if (BasisOptimize_print_lvl > 0) {
00721     writelogf("Error exit: returning with nComponents = %d\n", jr);
00722       }
00723 #endif
00724     return jr;
00725       }
00726      
00727       /*
00728        *  Assign a large negative number to the element that we have
00729        *  just found, in order to take it out of further consideration.
00730        */
00731       eAbund[kk] = test;
00732      
00733       /* *********************************************************** */
00734       /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX    */
00735       /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX  ****** */
00736       /* *********************************************************** */
00737       /*    
00738        *          Modified Gram-Schmidt Method, p. 202 Dalquist 
00739        *          QR factorization of a matrix without row pivoting. 
00740        */
00741       jl = jr;
00742       /*
00743        *   Fill in the row for the current element, k, under consideration
00744        *   The row will contain the Formula matrix value for that element
00745        *   with respect to the vector of component species.
00746        *   (note j and k indecises are flipped compared to the previous routine)
00747        */
00748       for (j = 0; j < nComponents; ++j) {
00749     jj = orderVectorSpecies[j];
00750     kk = orderVectorElements[k];
00751     sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
00752       } 
00753       if (jl > 0) {
00754     /*
00755      *         Compute the coefficients of JA column of the 
00756      *         the upper triangular R matrix, SS(J) = R_J_JR 
00757      *         (this is slightly different than Dalquist) 
00758      *         R_JA_JA = 1 
00759      */
00760     for (j = 0; j < jl; ++j) {
00761       ss[j] = 0.0;
00762       for (i = 0; i < nComponents; ++i) {
00763         ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
00764       }
00765       ss[j] /= sa[j];
00766     }
00767     /* 
00768      *     Now make the new column, (*,JR), orthogonal to the 
00769      *     previous columns
00770      */
00771     for (j = 0; j < jl; ++j) {
00772       for (l = 0; l < nComponents; ++l) {
00773         sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
00774       }
00775     }
00776       }
00777      
00778       /*
00779        *        Find the new length of the new column in Q. 
00780        *        It will be used in the denominator in future row calcs. 
00781        */
00782       sa[jr] = 0.0;
00783       for (ml = 0; ml < nComponents; ++ml) {
00784     double tmp = sm[ml + jr*nComponents];
00785     sa[jr] += tmp * tmp;
00786       }
00787       /* **************************************************** */
00788       /* **** IF NORM OF NEW ROW  .LT. 1E-6 REJECT ********** */
00789       /* **************************************************** */
00790       if (sa[jr] < 1.0e-6)  lindep = true;
00791       else                  lindep = false;
00792     } while(lindep);
00793     /* ****************************************** */
00794     /* **** REARRANGE THE DATA ****************** */
00795     /* ****************************************** */
00796     if (jr != k) {
00797 #ifdef DEBUG_MODE
00798       if (BasisOptimize_print_lvl > 0) {
00799         kk = orderVectorElements[k];
00800         ename = mphase->elementName(kk);
00801         writelog("   ---   ");
00802         writelogf("%-2.2s", ename.c_str()); 
00803         writelog("replaces ");
00804         kk = orderVectorElements[jr];
00805         ename = mphase->elementName(kk);
00806         writelogf("%-2.2s", ename.c_str()); 
00807         writelogf(" as element %3d\n", jr); 
00808       }
00809 #endif
00810       switch_pos(orderVectorElements, jr, k);
00811     }
00812   
00813     /*
00814      *      If we haven't found enough components, go back 
00815      *      and find some more. (nc -1 is used below, because
00816      *      jr is counted from 0, via the C convention.
00817      */
00818   } while (jr < (nComponents-1));
00819   return nComponents;
00820 } /* vcs_elem_rearrange() ****************************************************/
00821 
Generated by  doxygen 1.6.3