LinearEquation3< T > Class Template Reference

#include <LinearEquation3.h>

Classes
class	Array3

Public Member Functions
Basic3DVector< T >	solution (const Basic3DVector< T > &row0, const Basic3DVector< T > &row1, const Basic3DVector< T > &row2, const Basic3DVector< T > &rhsvec) const

Detailed Description

template<class T>
class LinearEquation3< T >

Definition at line 13 of file LinearEquation3.h.

Member Function Documentation

template<class T>

Basic3DVector<T> LinearEquation3< T >::solution ( const Basic3DVector< T > &  row0, const Basic3DVector< T > &  row1, const Basic3DVector< T > &  row2, const Basic3DVector< T > &  rhsvec  ) const

inline

Definition at line 55 of file LinearEquation3.h.

References funct::abs(), alignmentValidation::c1, gather_cfg::cout, FrontierConditions_GlobalTag_cff::dump, LinearEquation3< T >::Array3< U >::subtractScaled(), and std::swap().

Referenced by ThreePlaneCrossing::crossing().

                                                       {

    // copy the input to internal "matrix"
    Array3<T> row[3];
    row[0] = row0;
    row[1] = row1;
    row[2] = row2;
    Array3<T> rhs(rhsvec);

    // no implicit pivoting - rows expected to be normalized already

    // find pivot 0, i.e. row with largest first element
    int i0 = std::abs(row[0][0]) > std::abs(row[1][0]) ? 0 : 1;
    if (std::abs(row[i0][0]) < std::abs(row[2][0])) i0 = 2;

    int i1 = (i0+1)%3;
    int i2 = (i0+2)%3;

    // zero the first column of rows i1 and i2
    T c1 = row[i1][0] / row[i0][0];
    // row[i1] -= c1*row[i0];
    row[i1].subtractScaled( row[i0], c1);
    rhs[i1] -= c1*rhs[i0];
    T c2 = row[i2][0] / row[i0][0];
    // row[i2] -= c2*row[i0];
    row[i2].subtractScaled( row[i0], c2);
    rhs[i2] -= c2*rhs[i0];

    // find pivot 1, i.e. which row (i1 or i2) has the largest second element
    if (std::abs(row[i1][1]) < std::abs(row[i2][1])) std::swap( i1, i2);

    // zero the second column of row i2
    T c3 = row[i2][1] / row[i1][1];
    row[i2][1] -= c3 * row[i1][1];
    row[i2][2] -= c3 * row[i1][2];
    rhs[i2] -= c3*rhs[i1];

    // compute the solution
    T x2 = rhs[i2] / row[i2][2];
    T x1 = (rhs[i1] - x2*row[i1][2]) / row[i1][1];
    T x0 = (rhs[i0] - x1*row[i0][1] - x2*row[i0][2]) / row[i0][0];

    return Basic3DVector<T>(x0, x1, x2);
  }