умножить матрицы в петле в C

Question

Я пытаюсь умножить несколько матриц с помощью петли в C. Я получаю ожидаемый ответ в R, но не могу получить ожидаемый ответ в C. Я подозреваю, что проблема связана с функцией +=, которая, кажется, удваивает значение продукта после первой итерации цикла.

Я не очень хорошо знаком с C и не смог заменить функцию += той, которая вернет ожидаемый ответ.

Благодарим за любые советы.

Во-первых, здесь R код, который возвращает ожидаемый ответ:

B0 = -0.40 
B1 = 0.20 

mycov1 = exp(B0 + -2 * B1)/(1 + exp(B0 + -2 * B1)) 
mycov2 = exp(B0 + -1 * B1)/(1 + exp(B0 + -1 * B1)) 
mycov3 = exp(B0 + 0 * B1)/(1 + exp(B0 + 0 * B1)) 
mycov4 = exp(B0 + 1 * B1)/(1 + exp(B0 + 1 * B1)) 

trans1 = matrix(c(1 - 0.25 - mycov1,  mycov1, 0.25 * 0.80,    0, 
            0, 1 - 0.50,    0, 0.50 * 0.75, 
            0,   0,    1,    0, 
            0,   0,    0,    1), 
       nrow=4, ncol=4, byrow=TRUE) 

trans2 = matrix(c(1 - 0.25 - mycov2,  mycov2, 0.25 * 0.80,    0, 
            0, 1 - 0.50,    0, 0.50 * 0.75, 
            0,   0,    1,    0, 
            0,   0,    0,    1), 
       nrow=4, ncol=4, byrow=TRUE) 

trans3 = matrix(c(1 - 0.25 - mycov3,  mycov3, 0.25 * 0.80,    0, 
            0, 1 - 0.50,    0, 0.50 * 0.75, 
            0,   0,    1,    0, 
            0,   0,    0,    1), 
       nrow=4, ncol=4, byrow=TRUE) 

trans4 = matrix(c(1 - 0.25 - mycov4,  mycov4, 0.25 * 0.80,    0, 
            0, 1 - 0.50,    0, 0.50 * 0.75, 
            0,   0,    1,    0, 
            0,   0,    0,    1), 
       nrow=4, ncol=4, byrow=TRUE) 

trans2b <- trans1 %*% trans2 
trans3b <- trans2b %*% trans3 
trans4b <- trans3b %*% trans4 
trans4b 

# 
# This is the expected answer 
# 
#   [,1]  [,2]  [,3]  [,4] 
# [1,] 0.01819965 0.1399834 0.3349504 0.3173467 
# [2,] 0.00000000 0.0625000 0.0000000 0.7031250 
# [3,] 0.00000000 0.0000000 1.0000000 0.0000000 
# [4,] 0.00000000 0.0000000 0.0000000 1.0000000 
# 

Вот мой C код. C код довольно долго, потому что я не знаю C достаточно хорошо, чтобы быть эффективными:

#include <string.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 

char quit; 

int main(){ 

int i, j, k, ii, jj, kk ; 

double B0, B1, mycov ; 

double trans[4][4]  = {0} ; 

double prevtrans[4][4] = {{1,0,0,0}, 
          {0,1,0,0}, 
          {0,0,1,0}, 
          {0,0,0,1}}; 

B0 = -0.40 ; 
B1 = 0.20 ; 

for (i=1; i <= 4; i++) { 

      mycov = exp(B0 + B1 * (-2+i-1))/(1 + exp(B0 + B1 * (-2+i-1))) ; 

      trans[0][0] =  1 - 0.25 - mycov ; 
      trans[0][1] =    mycov ; 
      trans[0][2] =   0.25 * 0.80 ; 
      trans[0][3] =     0 ; 

      trans[1][0] =     0 ; 
      trans[1][1] =    1 - 0.50 ; 
      trans[1][2] =     0 ; 
      trans[1][3] =   0.50 * 0.75 ; 

      trans[2][0] =     0 ; 
      trans[2][1] =     0 ; 
      trans[2][2] =     1 ; 
      trans[2][3] =     0 ; 

      trans[3][0] =     0 ; 
      trans[3][1] =     0 ; 
      trans[3][2] =     0 ; 
      trans[3][3] =     1 ; 

      for (ii=0; ii<4; ii++){ 

       for(jj=0; jj<4; jj++){ 

        for(kk=0; kk<4; kk++){ 

         trans[ii][jj] += trans[ii][kk] * prevtrans[kk][jj] ; 

        } 
       } 
      } 

      prevtrans[0][0] =  trans[0][0] ; 
      prevtrans[0][1] =  trans[0][1] ; 
      prevtrans[0][2] =  trans[0][2] ; 
      prevtrans[0][3] =  trans[0][3] ; 
      prevtrans[1][0] =  trans[1][0] ; 
      prevtrans[1][1] =  trans[1][1] ; 
      prevtrans[1][2] =  trans[1][2] ; 
      prevtrans[1][3] =  trans[1][3] ; 
      prevtrans[2][0] =  trans[2][0] ; 
      prevtrans[2][1] =  trans[2][1] ; 
      prevtrans[2][2] =  trans[2][2] ; 
      prevtrans[2][3] =  trans[2][3] ; 
      prevtrans[3][0] =  trans[3][0] ; 
      prevtrans[3][1] =  trans[3][1] ; 
      prevtrans[3][2] =  trans[3][2] ; 
      prevtrans[3][3] =  trans[3][3] ; 

} 

    printf("To close this program type 'quit' and hit the return key\n"); 
    printf("    \n"); 
    scanf("%d", &quit); 

    return 0; 

} 

Вот окончательная матрица возвращаемый выше C код:

0.4821 3.5870 11.68 381.22 
0  1  0  76.875 
0  0  5  0 
0  0  0  5 

Answer 1

Эта линия

     trans[ii][jj] += trans[ii][kk] * prevtrans[kk][jj] ; 

не в порядке. Вы изменяете trans на месте, пока вы все еще используете его для вычисления результирующей матрицы. Вам понадобится другая матрица для временного хранения результата умножения. И затем используйте:

 // Store the resultant matrix in temp. 
     for (ii=0; ii<4; ii++){ 

      for(jj=0; jj<4; jj++){ 

       temp[ii][jj] = 0.0; 

       for(kk=0; kk<4; kk++){ 

        temp[ii][jj] += trans[ii][kk] * prevtrans[kk][jj] ; 

       } 
      } 
     } 

     // Transfer the data from temp to trans 
     for (ii=0; ii<4; ii++){ 

      for(jj=0; jj<4; jj++){ 
       trans[ii][jj] = temp[ii][jj]; 
      } 
     }