Get the right data type from MatrixExp?

Questions
1

Good morning,
I am using MatrixExp[] inside a Compiled function and, at a certain point, I assign to a square matrix the output of MatrixExp. I am getting an error complaining that MatrixExp is of type {_Real,0} and cannot be assigned to a variable of type {_Real,2}. The code is quite straightforward… The problem occurs with the instruction “M=M.Mexp”. If I change that instruction to “M = M + Mexp” the problem goes away. Is M.Mexp not of type {Real,2}?? Any thoughts?

par = {5.7, 1.038, 1.016, 1.018, 1.015, 0.186, 0.003, 0.009, 0.004, 
   0.0003, 0.00004, 0.071};
dim = 69;
P = Compile[{{par, _Real, 1}},
   Module[{i, L1, L2, L3, L4, L, q1, q2, q3, q4, gama, alfa, beta, 
     A},
    L1 = par[[1]]; L2 = par[[2]]; L3 = par[[3]]; L4 = par[[4]]; 
    L = par[[5]];
    q1 = par[[6]]; q2 = par[[7]]; q3 = par[[8]]; q4 = par[[9]]; 
    gama = par[[10]];
    alfa = par[[11]]; beta = par[[12]];
    A = ConstantArray[0.0, {dim, dim}];
    A[[1, 1]] = -(L1 + q1); A[[2, 2]] = -(L2 + q2); 
    A[[3, 3]] = -(L3 + q3); A[[4, 4]] = -(L4 + q4);
    A[[1, 2]] = L1; A[[2, 3]] = L2; A[[3, 4]] = L3; A[[4, 5]] = L4;
    For[i = 5, i <= dim - 1, i++, A[[i, i + 1]] = L];
    For[i = 5, i <= dim - 1, i++, 
     A[[i, i]] = -(L + gama + alfa*Exp[beta*i])];
    A[[dim, dim]] = -(gama + alfa*Exp[beta*dim]);
    Return[A]
    ]
   ];
OBJ2 = Compile[{{par, _Real, 1}},
   Module[{M, Mexp, s, i},
    Mexp = MatrixExp[P[par]];
    M = ConstantArray[0., {dim, dim}];
    For[i = 1, i <= dim, i++, M[[i, i]] = 1.;];
    vec = ConstantArray[0., {dim + 1}];
    For[i = 1, i <= dim, i++,
     M = M.Mexp; vec[[i]] = Total[M[[1, All]]]];
    
    Return[
     Total[Table[((Q[[i]] - (1 - vec[[i + 1]]/vec[[i]]))^2)*
        S[[i]]^2, {i, 1, dim}]]]
    ]
   ];
2

As it stands it has also errors with the Q[] and S[]

Clear[OBJ2]
OBJ2 = Compile[{{par, _Real, 1}},
  Module[{M, Mexp, s, i},
   Mexp = MatrixExp[P[par]];
   M = ConstantArray[0., {dim, dim}];
   For[i = 1, i <= dim, i++,
    M[[i, i]] = 1.
    ];
   vec = ConstantArray[0., {dim + 1}];
   For[i = 1, i <= dim, i++,
    M = M.Mexp;
    vec[[i]] = Total[M[[1, All]]]
    ];
   Return[
    Total[Table[((Q[[i]] - (1 - vec[[i + 1]]/vec[[i]]))^2)*
       S[[i]]^2, {i, 1, dim}]]]
   ]
  ]

Compile::cset: Variable M of type {_Real,2} encountered in assignment of type {_Real,0}.

Compile::part: Part specification Q[[i]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

Compile::part: Part specification S[[i]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

Compile::part: Part specification Q[[i]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

General::stop: Further output of Compile::part will be suppressed during this calculation.

Compile::cset: Variable M of type {_Real,2} encountered in assignment of type {_Real,0}.

I skipped the S[]and Q[], gave the functions all arguments (dim), all variables (vec) and cleared an useless for loop by using IdentityMatrix[]. Still the error with the tensor rank appears under Compile[] and is of course a wrong error message. You might file a bug to Wolfram.

But the good message is, with the function o1 you can compute your stuff rather fast without compile (from o1 we know that Compile[] has a bug). See the notebook attached for details.