to perform this += T' * B * T
| 943 | |
| 944 | // to perform this += T' * B * T |
| 945 | int |
| 946 | Matrix::addMatrixTripleProduct(double thisFact, |
| 947 | const Matrix &T, |
| 948 | const Matrix &B, |
| 949 | double otherFact) |
| 950 | { |
| 951 | if (thisFact == 1.0 && otherFact == 0.0) |
| 952 | return 0; |
| 953 | #ifdef _G3DEBUG |
| 954 | if ((numCols != numRows) || (B.numCols != B.numRows) || (T.numCols != numRows) || |
| 955 | (T.numRows != B.numCols)) { |
| 956 | opserr << "Matrix::addMatrixTripleProduct() - incompatable matrices\n"; |
| 957 | return -1; |
| 958 | } |
| 959 | #endif |
| 960 | |
| 961 | // cheack work area can hold the temporary matrix |
| 962 | int dimB = B.numCols; |
| 963 | int sizeWork = dimB * numCols; |
| 964 | |
| 965 | if (sizeWork > sizeDoubleWork) { |
| 966 | this->addMatrix(thisFact, T^B*T, otherFact); |
| 967 | return 0; |
| 968 | } |
| 969 | |
| 970 | // zero out the work area |
| 971 | double *matrixWorkPtr = matrixWork; |
| 972 | for (int l=0; l<sizeWork; l++) |
| 973 | *matrixWorkPtr++ = 0.0; |
| 974 | |
| 975 | // now form B * T * fact store in matrixWork == A area |
| 976 | // NOTE: looping as per blas3 dgemm_: j,k,i |
| 977 | |
| 978 | double *tkjPtr = &(T.data)[0]; |
| 979 | for (int j=0; j<numCols; j++) { |
| 980 | double *aijPtrA = &matrixWork[j*dimB]; |
| 981 | for (int k=0; k<dimB; k++) { |
| 982 | double tmp = *tkjPtr++ * otherFact; |
| 983 | double *aijPtr = aijPtrA; |
| 984 | double *bikPtr = &(B.data)[k*dimB]; |
| 985 | for (int i=0; i<dimB; i++) |
| 986 | *aijPtr++ += *bikPtr++ * tmp; |
| 987 | } |
| 988 | } |
| 989 | |
| 990 | // now form T' * matrixWork |
| 991 | // NOTE: looping as per blas3 dgemm_: j,i,k |
| 992 | if (thisFact == 1.0) { |
| 993 | double *dataPtr = &data[0]; |
| 994 | for (int j=0; j< numCols; j++) { |
| 995 | double *workkjPtrA = &matrixWork[j*dimB]; |
| 996 | for (int i=0; i<numRows; i++) { |
| 997 | double *ckiPtr = &(T.data)[i*dimB]; |
| 998 | double *workkjPtr = workkjPtrA; |
| 999 | double aij = 0.0; |
| 1000 | for (int k=0; k< dimB; k++) |
| 1001 | aij += *ckiPtr++ * *workkjPtr++; |
| 1002 | *dataPtr++ += aij; |