LAPACK QR-факторизация

Я пытаюсь получить Q-Matrix из подпрограмм LAPACK zgeqrf и zungqr.

У меня есть комплексная матрица Nw-by-Nw с неортогональными векторами в ее столбцах. Это мой код на C ++. (Матрица называется vr_tr)

//QR Fact.
complex<double> TAU[Nw*Nw];
zgeqrf_(&Nw, &Nw, vr_tr, &Nw, TAU, &wkopt, &lwork, &info);
lwork = (int)wkopt.real();
work = new complex<double> [lwork];
zgeqrf_(&Nw, &Nw, vr_tr, &Nw, TAU, work, &lwork, &info);
zungqr_(&Nw, &Nw, &Nw, vr_tr, &Nw, TAU, &wkopt, &lwork, &info);
lwork = (int)wkopt.real();
work = new complex<double> [lwork];
zungqr_(&Nw, &Nw, &Nw, vr_tr, &Nw, TAU, work, &lwork, &info);

//Checking if vr_tr * vr_tr' = I
complex<double> one = 1,zer=0,res[Nw*Nw]; char Tchar = 'T', Nchar = 'N';
zgemm_( &Nchar, &Tchar, &Nw, &Nw, &Nw, &one, vr_tr, &Nw, vr_tr, &Nw, &zer, res, &Nw );
for(i=0;i<Nw;i++){
    for(j=0;j<Nw;j++){
        cout<<res[i*Nw+j]<<"\t";
    }
    cout<<"\n\n";
}

То, что я получаю после запуска этого кода, не является единичной матрицей, как предполагалось, потому что он должен вычислять QR-факторизацию из zgeqrf и получать Q-матрицу из zungqr, а Q-матрица должна быть ортогональной, поэтому Q*Q'=I.

Что не так с этим кодом?


c++ matrix lapack
person Alireza    schedule 07.08.2017    source источник
comment
Из любопытства, зачем вы возитесь с лапаком, когда существуют высокоуровневые оболочки C ++, такие как Armadillo ?   -  person The Quantum Physicist    schedule 07.08.2017
comment
@TheQuantumPhysicist Habits и потому, что люди, с которыми я работаю, используют его. Однако никогда не пробовал Армадилло. Может ли он выполнять все процедуры LAPACK?   -  person Alireza    schedule 07.08.2017
comment
Проверьте документацию по предоставленной мной ссылке. Когда я все еще занимался исследованиями в академических кругах, там были все необходимые мне процедуры и даже больше. Например, LAPACK не имеет универсальной функции возведения в степень матрицы (если вы не хотите делать это самостоятельно, вычисляя собственные значения и предполагая, что у вас есть симметричная матрица), но у Armadillo она есть. Вы можете связать Armadillo с любой реализацией LAPACK по вашему желанию, и это сэкономит вам много времени. Единственная проблема с Armadillo заключается в том, что он не поддерживает кластеры (AFAIK). В противном случае вам следует пойти на это.   -  person The Quantum Physicist    schedule 07.08.2017
comment
@TheQuantumPhysicist Отлично! Спасибо за совет!   -  person Alireza    schedule 07.08.2017

Ответы (1)


arrow_upward
0
arrow_downward

Я использовал zunmqr вместо znugqr, и это было исправлено.

Хотя, честно говоря, не знаю, почему zungqr не работал.

person Alireza    schedule 07.08.2017