Method addMatrixTransposeVector

SRC/matrix/Vector.cpp:459–599

int
Vector::addMatrixTransposeVector(double thisFact, 
				 const Matrix &m, 
				 const Vector &v, 
				 double otherFact )
{
  // see if quick return
  if (otherFact == 0.0 && thisFact == 1.0)
    return 0;

#ifdef _G3DEBUG
  // check the sizes are compatible
  if ((sz != m.noRows()) && (m.noRows() != v.sz)) {
    // otherwise incompatible sizes
    opserr << "Vector::addMatrixTransposeVector() - incompatible sizes\n";
    return -1;    
  }
#endif

  if (thisFact == 1.0) {

    // want: this += m^t * v * otherFact
    if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
      int otherSize = v.sz;
      double *matrixDataPtr = m.data;
      double *otherDataPtrA = v.theData;
      for (int i=0; i<sz; i++) {
	double *otherDataPtr = otherDataPtrA;
	double sum = 0.0;
	for (int j=0; j<otherSize; j++)
	  sum += *matrixDataPtr++ * *otherDataPtr++;
	theData[i] += sum;
      }
    } else if (otherFact == -1.0) { // no point doing multiplication if otherFact = 1.0
      int otherSize = v.sz;
      double *matrixDataPtr = m.data;
      double *otherDataPtrA = v.theData;
      for (int i=0; i<sz; i++) {
	double *otherDataPtr = otherDataPtrA;
	double sum = 0.0;
	for (int j=0; j<otherSize; j++)
	  sum += *matrixDataPtr++ * *otherDataPtr++;
	theData[i] -= sum;
      }
    } else {
      int otherSize = v.sz;
      double *matrixDataPtr = m.data;
      double *otherDataPtrA = v.theData;
      for (int i=0; i<sz; i++) {
	double *otherDataPtr = otherDataPtrA;
	double sum = 0.0;
	for (int j=0; j<otherSize; j++)
	  sum += *matrixDataPtr++ * *otherDataPtr++;
	theData[i] += sum * otherFact;
      }
    }
  }