[РЕШЕНО] ALGLIB и многопоточность

[Vamireh]

Решил избавиться от своего класса решения СЛАУ и использовать ALGLIB. СЛАУ у меня решается тысячами, поэтому многопоточно (через concurrent). В итоге получаю, что первая система решается нормально, а все остальные - там фиг знает что. Если решать в один поток в цикле, то все нормально. Кто-нибудь сталкивался с чем-то подобным?

[Vamireh]

Да, еще. Если запускать с дебаггером, то все хорошо.

[Old]

Решил избавиться от своего класса решения СЛАУ и использовать ALGLIB. СЛАУ у меня решается тысячами, поэтому многопоточно (через concurrent). В итоге получаю, что первая система решается нормально, а все остальные - там фиг знает что. Если решать в один поток в цикле, то все нормально. Кто-нибудь сталкивался с чем-то подобным?

Так а где код?

[Vamireh]

Код:

 Kramer3d kramer;
    auto end = params.data->constEnd();

    real_2d_array a; a.setlength(3, 3);
    real_1d_array b; b.setlength(3);

    kramer.a[0][0] = Size;
    a(0,0) = Size;
    for (auto it = params.data->constBegin(); it != end; ++it) {
        const double tdTau = it->t / params.tau;
        kramer.b[0] += it->NR[params.axe];
        kramer.b[1] += it->t * it->NR[params.axe];
        kramer.b[2] += qExp(-tdTau) * it->NR[params.axe];
        kramer.a[0][1] += it->t;
        kramer.a[0][2] += qExp(-tdTau);
        kramer.a[1][1] += qPow(it->t, 2.0);
        kramer.a[1][2] += qExp(-tdTau) * it->t;
        kramer.a[2][2] += qExp(-2.0 * tdTau);

        b(0) += it->NR[params.axe];
        b(1) += it->t * it->NR[params.axe];
        b(2) += qExp(-tdTau) * it->NR[params.axe];
        a(0, 1) += it->t;
        a(0, 2) += qExp(-tdTau);
        a(1, 1) += qPow(it->t, 2.0);
        a(1, 2) += qExp(-tdTau) * it->t;
        a(2, 2) += qExp(-2.0 * tdTau);
    }
    kramer.a[1][0] = kramer.a[0][1];
    kramer.a[2][0] = kramer.a[0][2];
    kramer.a[2][1] = kramer.a[1][2];

    a(1, 0) = a(0, 1);
    a(2, 0) = a(0, 2);
    a(2, 1) = a(1, 2);

    const Array3D slu = kramer.calcSLU();
    result.koefs = {slu[2], slu[1], slu[0], static_cast<double>(params.tau)};
    real_1d_array root;
    ae_int_t info;
    densesolverreport rep;
    rmatrixsolve(a, 3, b, info, rep, root);

    if (params.axe == 0 && params.tau < 100)
    qDebug() << root(0) << root(1) << root(2) << "\t-\t" << slu[0] << slu[1] << slu[2];

kramer - это мое.

[Old]

Как запускается многопоточный расчет?

[Vamireh]

Как запускается многопоточный расчет?

Код:

Vector<OneAxeParams> paramsOneAxe;
        int tau = Devices.at(params.deviceTypeIndex).tauMin;
        while (tau <= Devices.at(params.deviceTypeIndex).tauMax)
            paramsOneAxe.append( {&internalData.calcedData, axe, tau++} );

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
        internalData.resultsForAllTau[axe] = QtConcurrent::blockingMapped(paramsOneAxe, &RealizInfo::calcOneAxe);
#pragma GCC diagnostic pop

[Vamireh]

еще заметил, что в ALGLIB в a элементы матрицы текущей итерации иногда равны сумме с предыдущего шага и с текущего шага, а иногда они сбрасываются и снова потом прибавляются. Имеются ввиду правильные значения, из моего класса.

[Old]

Можно еще посмотреть что из себя представляет internalData, OneAxeParams и RealizInfo::calcOneAxe

И что за warning вы отключили? А то не понятно, что где затеняет?

[Fat-Zer]

http://www.alglib.net/
как я понял, gpl'ная версия однопоточная... угадал?

[Vamireh]

В метро ехал, меня осенило. Я - то в конструкторе Kramer'а обнуляю поля, а ALGLIB - нет, завтра проверю. Но сейчас дома примитивнейший пример аналогичный написал - все работает.

[Old]

В метро ехал, меня осенило. Я - то в конструкторе Kramer'а обнуляю поля, а ALGLIB - нет, завтра проверю. Но сейчас дома примитивнейший пример аналогичный написал - все работает.

Когда в debug работает, а в release нет - не инициализированные данные причина номер один.

[Vamireh]

Это да, просто я в своих классах всегда инициирую поля, а тут по привычке не подумал сразу об этом.