[kernel] direct fill HTrans with the whole block
@@ -531,6 +531,8 @@ void OSNSMatrix::fillH(DynamicalSystemsGraph & DSG, InteractionsGraph& indexSet,
}
DEBUG_END("void OSNSMatrix::fillH(SP::DynamicalSystemsGraph DSG, InteractionsGraph& indexSet, bool update)\n");
}
+
+#include <float.h>
// Fill the matrix H
void OSNSMatrix::fillHtrans(DynamicalSystemsGraph & DSG, InteractionsGraph& indexSet, bool update)
{
@@ -555,32 +557,34 @@ void OSNSMatrix::fillHtrans(DynamicalSystemsGraph & DSG, InteractionsGraph& inde
NM_triplet_alloc(&H_NM, _triplet_nzmax);
CSparseMatrix* Htriplet= NM_triplet(&H_NM);
-
- unsigned int pos = 0, pos_ds=0;
+ unsigned int pos = 0, abs_pos_ds=0;
SP::SiconosMatrix leftInteractionBlock;
+
+
InteractionsGraph::VIterator ui, uiend;
for(std::tie(ui, uiend) = indexSet.vertices(); ui != uiend; ++ui)
{
Interaction& inter = *indexSet.bundle(*ui);
+ size_t sizeY = inter.dimension();
+ leftInteractionBlock = inter.getLeftInteractionBlock();
+
+ double * array = &*leftInteractionBlock->getArray();
+
SP::DynamicalSystem ds1 = indexSet.properties(*ui).source;
SP::DynamicalSystem ds2 = indexSet.properties(*ui).target;
bool endl = false;
size_t posBlock = indexSet.properties(*ui).source_pos;
- size_t pos2 = indexSet.properties(*ui).target_pos;
+ size_t pos_ds2 = indexSet.properties(*ui).target_pos;
pos = indexSet.properties(*ui).absolute_position;
- for(SP::DynamicalSystem ds = ds1; !endl; ds = ds2, posBlock = pos2)
+
+ for(SP::DynamicalSystem ds = ds1; !endl; ds = ds2, posBlock = pos_ds2)
{
endl = (ds == ds2);
size_t sizeDS = ds->dimension();
- size_t sizeY = inter.dimension();
- // this whole part is a hack. Just should just get the rightblock
- leftInteractionBlock = inter.getLeftInteractionBlockForDS(posBlock, sizeY, sizeDS);
- pos_ds = DSG.properties(DSG.descriptor(ds)).absolute_position;
- DEBUG_PRINTF("pos = %u", pos);
- DEBUG_PRINTF("pos_ds = %u", pos_ds);
- leftInteractionBlock->fillTriplet(Htriplet, pos, pos_ds);
+ abs_pos_ds = DSG.properties(DSG.descriptor(ds)).absolute_position;
+ CSparseMatrix_block_dense_zentry(Htriplet, pos, abs_pos_ds, array+posBlock*sizeY, sizeY, sizeDS, DBL_EPSILON);
}
}
_triplet_nzmax = NM_nnz(&H_NM);
