Наименьшие квадраты с использованием scalaPACK ⇐ C++

[Anonymous]

Я учусь использовать процедуру scalaPACK для метода наименьших квадратов, но не могу заставить ее работать.
Буду признателен за ваше понимание этого вопроса.
В коде ниже предполагается распределение матрицы A 8x8 на 4 локальные матрицы A_loc 4x4 и аналогично для вектора B. Черная сетка имеет размер 2x2. Код предназначен для запуска с четырьмя задачами MPI.
Код компилируется и запускается со следующим выводом с одной задачей MPI:
BEFORE: Local B size 8 on rank: 0 values: 1 2 3 4 5 6 7 8
AFTER: Local B size 8 on rank: 0 values: 0.960918 0.0252256 0.0222086 0.0319953 -0.0471391 0.0205009 -0.0163878 0.00267855

Выглядит правильно, B_loc был обновлен решением по методу наименьших квадратов.
Однако, когда используются четыре задачи MPI, я получаю следующий результат:
BEFORE: Local B size 4 on rank: 0 values: 1 2 3 4
BEFORE: Local B size 4 on rank: 1 values: 4.68603e-310 0 2.122e-314 2.122e-314
BEFORE: Local B size 4 on rank: 2 values: 5 6 7 8
BEFORE: Local B size 4 on rank: 3 values: 4.67566e-310 0 2.122e-314 2.122e-314
AFTER: Local B size 4 on rank: 1 values: 4.68603e-310 0 2.122e-314 2.122e-314
AFTER: Local B size 4 on rank: 2 values: 5 6 7 8
AFTER: Local B size 4 on rank: 3 values: 4.67566e-310 0 2.122e-314 2.122e-314
AFTER: Local B size 4 on rank: 0 values: 0.321957 1.12458 -0.21503 -0.231507

Что меня удивляет, так это то, что только значения B_loc в ранге 0 были обновлены, предположительно, с использованием части решения методом наименьших квадратов. Что я делаю не так?
Вот код:
#include
#include
#include

extern "C" void blacs_get_(int*, int*, int*);
extern "C" void blacs_pinfo_(int*, int*);
extern "C" void blacs_gridinit_(int*, char*, int*, int*);
extern "C" void blacs_gridinfo_(int*, int*, int*, int*, int*);
extern "C" void descinit_(int*, int*, int*, int*, int*, int*, int*, int*, int*, int*);
extern "C" void blacs_gridexit_(int*);
extern "C" int numroc_(int*, int*, int*, int*, int*);
extern "C" void pdgels_(char* trans, int* m, int* n, int* nrhs, double* a, int* ia, int* ja, int* desca, double* b, int* ib, int* jb, int* descb, double* work, int* lwork, int* info);
extern "C" void pdgemr2d_(int *m, int *n, double *a, int *ia, int *ja, int *desca,
double *b, int *ib, int *jb, int *descb, int *context);

int main() {

MPI_Init(NULL, NULL);
int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);

int b_nrows, b_ncols, b_row, b_col, myid;
blacs_pinfo_( &myid, &size );
char layout = 'R';

b_nrows = sqrt(size);
b_ncols = sqrt(size);

int context;
int izero=0;
int ione=1;
blacs_get_( &izero, &izero, &context );
blacs_gridinit_( &context, &layout, &b_nrows, &b_ncols );
blacs_gridinfo_( &context, &b_nrows, &b_ncols, &b_row, &b_col );

int m=8; // rows of the global array
int n=8; // cols
int mb=m;
int nb=n;
int mb_loc=m/b_nrows;
int nb_loc=n/b_ncols;

// size of local arrays
int m_loc = numroc_( &m, &mb_loc, &b_row, &izero, &b_nrows );
int n_loc = numroc_( &n, &nb_loc, &b_col, &izero, &b_ncols );

// global
double *A = new double[m*n] ;
double *B = new double[m] ;
// local
double *A_loc = new double[m_loc*n_loc] ;
double *B_loc = new double[m_loc] ;

int descA[9], descB[9], descA_loc[9], descB_loc[9];
int info;
descinit_( descA, &m, &n, &mb, &nb, &izero, &izero, &context, &m, &info);
descinit_( descB, &m, &ione, &mb, &ione, &izero, &izero, &context, &m, &info);

descinit_( descA_loc, &m, &n, &mb_loc, &nb_loc, &izero, &izero, &context, &m_loc, &info);
descinit_( descB_loc, &m, &ione, &mb_loc, &ione, &izero, &izero, &context, &m_loc, &info);

// generate column major order matrix
for (int i = 0; i < m*n; i++) {
A = i+1;
}
// and vector of tagets
for (int i = 0; i < m; i++) {
B = (i+1);
}

// distribute matrices
pdgemr2d_(&m, &n, A, &ione, &ione, descA,
A_loc, &ione, &ione, descA_loc, &context);

pdgemr2d_(&m, &ione, B, &ione, &ione, descB,
B_loc, &ione, &ione, descB_loc, &context);

std::cout

Подробнее здесь: https://stackoverflow.com/questions/783 ... -scalapack