#include "var1EXX.h"

using namespace scl;
using namespace std;

void var1EXX
(const realmat& b0, const realmat& B, const realmat& S, realmat& EXX)
{
   INTEGER M = b0.size();
   if (B.nrow() != M || B.ncol() != M) error("Error, wrong dimensions");
   if (S.nrow() != M || S.ncol() != M) error("Error, wrong dimensions");

   realmat I(M,M,0.0);
   for (INTEGER i=1; i<=M; ++i) I(i,i) = 1.0;

   realmat IB = I - B;
   realmat IBinv = inv(IB);

   realmat mu = IBinv*b0;

   realmat vecS(M*M,1);
   for (INTEGER i=1; i<=M*M; ++i) vecS[i] = S[i];

   realmat vecV = inv(kronprd(I,I) - kronprd(B,B))*vecS;

   realmat C22(M,M);

   for (INTEGER i=1; i<=M*M; ++i) C22[i] = vecV[i];

   for (INTEGER j=1; j<=M; ++j) {
     for (INTEGER i=1; i<=M; ++i) {
       C22(i,j) += mu[i]*mu[j];
     }
   }

   EXX.resize(M+1,M+1);

   EXX(1,1) = 1.0;

   for (INTEGER i=1; i<=M; ++i) EXX(1+i,1) = EXX(1,1+i) = mu[i];

   for (INTEGER j=1; j<=M; ++j) {
     for (INTEGER i=1; i<=M; ++i) {
       EXX(1+i,1+j) = C22(i,j);
     }
   }
}
