Method addMatrixTripleProduct

SRC/matrix/Matrix.cpp:936–1027

to perform this += T' * B * T

// to perform this += T' * B * T
int
Matrix::addMatrixTripleProduct(double thisFact, 
			       const Matrix &T, 
			       const Matrix &B, 
			       double otherFact)
{
    if (thisFact == 1.0 && otherFact == 0.0)
      return 0;
#ifdef _G3DEBUG
    if ((numCols != numRows) || (B.numCols != B.numRows) || (T.numCols != numRows) ||
	(T.numRows != B.numCols)) {
      opserr << "Matrix::addMatrixTripleProduct() - incompatible matrices\n";
      return -1;
    }
#endif

    // cheack work area can hold the temporary matrix
    int dimB = B.numCols;
    int sizeWork = dimB * numCols;

    if (sizeWork > sizeDoubleWork) {
      this->addMatrix(thisFact, T^B*T, otherFact);
      return 0;
    }

    // zero out the work area
    double *matrixWorkPtr = matrixWork;
    for (int l=0; l<sizeWork; l++)
      *matrixWorkPtr++ = 0.0;
    
    // now form B * T * fact store in matrixWork == A area
    // NOTE: looping as per blas3 dgemm_: j,k,i

    double *tkjPtr  = &(T.data)[0];
    for (int j=0; j<numCols; j++) {
      double *aijPtrA = &matrixWork[j*dimB];
      for (int k=0; k<dimB; k++) {
	double tmp = *tkjPtr++ * otherFact;
	double *aijPtr = aijPtrA;
	double *bikPtr = &(B.data)[k*dimB];
	for (int i=0; i<dimB; i++) 
	  *aijPtr++ += *bikPtr++ * tmp;
      }
    }

    // now form T' * matrixWork
    // NOTE: looping as per blas3 dgemm_: j,i,k
    if (thisFact == 1.0) {
      double *dataPtr = &data[0];
      for (int j=0; j< numCols; j++) {
	double *workkjPtrA = &matrixWork[j*dimB];
	for (int i=0; i<numRows; i++) {
	  double *ckiPtr = &(T.data)[i*dimB];
	  double *workkjPtr = workkjPtrA;
	  double aij = 0.0;
	  for (int k=0; k< dimB; k++)
	    aij += *ckiPtr++ * *workkjPtr++;
	  *dataPtr++ += aij;