Commit 28f2f1f3 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

ignore zeros entries in matrix/matrix multiplications in dirdynared

parent fb79f84e
......@@ -76,19 +76,27 @@ int dirdynared(MbsAux *mbs_aux, MbsData *s)
// Fr_uc = Fuc + (Muc_v+ Bvuc'*Mvv)*bprim + Bvuc'*Fv ;
// BtMvu = Bvuc'*Mv_uc
nL = nC = mbs_aux->iquc[0];
nk = s->nqv;
for (i = 1; i <= nL; i++)
{
for (j = 1; j <= nC; j++)
for (j = 1; j <= nC; j++)
mbs_aux->BtMvu[i][j] = 0;
}
for (k = 1; k <= nk; k++)
{
for (j = 1; j <= nC; j++)
{
double m = mbs_aux->M[s->qv[k]][mbs_aux->iquc[j]];
if (m == 0) continue;
for (i = 1; i <= nL; i++)
{
term = 0.0;
for (k = 1; k <= nk; k++)
{
term += mbs_aux->Bvuc[k][i] * mbs_aux->M[s->qv[k]][mbs_aux->iquc[j]];
}
mbs_aux->BtMvu[i][j] = term;
mbs_aux->BtMvu[i][j] += mbs_aux->Bvuc[k][i] *m;
}
}
}
// BtMvv = Bvuc'*Mvv
......@@ -96,15 +104,20 @@ int dirdynared(MbsAux *mbs_aux, MbsData *s)
nC = nk = s->nqv;
for (i = 1; i <= nL; i++)
{
for (j = 1; j <= nC; j++)
for (j = 1; j <= nC; j++)
mbs_aux->BtMvv[i][j] = 0;
}
for (j = 1; j <= nC; j++)
{
for (k = 1; k <= nk; k++)
{
double m = mbs_aux->M[s->qv[k]][s->qv[j]];
if (m == 0) continue;
for (i = 1; i <= nL; i++)
{
term = 0.0;
for (k = 1; k <= nk; k++)
{
term += mbs_aux->Bvuc[k][i] * mbs_aux->M[s->qv[k]][s->qv[j]];
}
mbs_aux->BtMvv[i][j] = term;
mbs_aux->BtMvv[i][j] += mbs_aux->Bvuc[k][i] * mbs_aux->M[s->qv[k]][s->qv[j]];
}
}
}
// BtFv = Bvuc'*Fv
......@@ -125,15 +138,20 @@ int dirdynared(MbsAux *mbs_aux, MbsData *s)
nk = s->nqv;
for (i = 1; i <= nL; i++)
{
for (j = 1; j <= nC; j++)
mbs_aux->BtMB[i][j] = 0;
}
for (i = 1; i <= nL; i++)
{
for (k = 1; k <= nk; k++)
{
double bm = mbs_aux->BtMvv[i][k];
if (bm == 0) continue;
for (j = 1; j <= nC; j++)
{
term = 0.0;
for (k = 1; k <= nk; k++)
{
term += mbs_aux->BtMvv[i][k] * mbs_aux->Bvuc[k][j];
}
mbs_aux->BtMB[i][j] = term;
mbs_aux->BtMB[i][j] += bm * mbs_aux->Bvuc[k][j];
}
}
}
// MBMb = (Muv+Bvuc'*Mvv)*bprim
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment