| 794 | |
| 795 | |
| 796 | int |
| 797 | Matrix::addMatrixProduct(double thisFact, |
| 798 | const Matrix &B, |
| 799 | const Matrix &C, |
| 800 | double otherFact) |
| 801 | { |
| 802 | if (thisFact == 1.0 && otherFact == 0.0) |
| 803 | return 0; |
| 804 | #ifdef _G3DEBUG |
| 805 | if ((B.numRows != numRows) || (C.numCols != numCols) || (B.numCols != C.numRows)) { |
| 806 | opserr << "Matrix::addMatrixProduct(): incompatible matrices, this\n"; |
| 807 | return -1; |
| 808 | } |
| 809 | #endif |
| 810 | // NOTE: looping as per blas3 dgemm_: j,k,i |
| 811 | if (thisFact == 1.0) { |
| 812 | |
| 813 | // want: this += B * C otherFact |
| 814 | int numColB = B.numCols; |
| 815 | double *ckjPtr = &(C.data)[0]; |
| 816 | for (int j=0; j<numCols; j++) { |
| 817 | double *aijPtrA = &data[j*numRows]; |
| 818 | for (int k=0; k<numColB; k++) { |
| 819 | double tmp = *ckjPtr++ * otherFact; |
| 820 | double *aijPtr = aijPtrA; |
| 821 | double *bikPtr = &(B.data)[k*numRows]; |
| 822 | for (int i=0; i<numRows; i++) |
| 823 | *aijPtr++ += *bikPtr++ * tmp; |
| 824 | } |
| 825 | } |
| 826 | } |
| 827 | |
| 828 | else if (thisFact == 0.0) { |
| 829 | |
| 830 | // want: this = B * C otherFact |
| 831 | double *dataPtr = data; |
| 832 | for (int i=0; i<dataSize; i++) |
| 833 | *dataPtr++ = 0.0; |
| 834 | int numColB = B.numCols; |
| 835 | double *ckjPtr = &(C.data)[0]; |
| 836 | for (int j=0; j<numCols; j++) { |
| 837 | double *aijPtrA = &data[j*numRows]; |
| 838 | for (int k=0; k<numColB; k++) { |
| 839 | double tmp = *ckjPtr++ * otherFact; |
| 840 | double *aijPtr = aijPtrA; |
| 841 | double *bikPtr = &(B.data)[k*numRows]; |
| 842 | for (int i=0; i<numRows; i++) |
| 843 | *aijPtr++ += *bikPtr++ * tmp; |
| 844 | } |
| 845 | } |
| 846 | } |
| 847 | |
| 848 | else { |
| 849 | // want: this = B * C otherFact |
| 850 | double *dataPtr = data; |
| 851 | for (int i=0; i<dataSize; i++) |
| 852 | *dataPtr++ *= thisFact; |
| 853 | int numColB = B.numCols; |
no outgoing calls
no test coverage detected