MKL ScaLapack problems

MKL ScaLapack problems

本文关键字:problems ScaLapack MKL      更新时间:2023-10-16

我试图从http://acts.nersc.gov/scalapack/hands-on/etc/pddttrdrv/pddttrdrv.c.html运行一个简单的Hello World (MKL) ScaLapack示例,但我偶然发现了一个问题(我使用MPICH2和我的操作系统是Windows)。

当我运行带有MPI标志的代码时

-localonly 2

:

{   -1,   -1}:  On entry to
{   -1,   -1}:  On entry to
 parameter number    1 had an illegal value
 parameter number    1 had an illegal value

PDDTTRF, D&C alg.: only 1 block per proc

我仔细检查了参数,并与在线参考文献进行了比较,看看这些值是否正确,我似乎没有发现问题。

代码是:

#include <mpi.h>
#define numroc_ NUMROC
#define descinit_ DESCINIT
#include <iostream>
#include <math.h>
#include <mkl_pblas.h>
#include <mkl_scalapack.h>
#include <mkl_blacs.h>
using namespace std;
extern "C" {
    /* Cblacs declarations */
    void Cblacs_pinfo(int*, int*);
    void Cblacs_get(int, int, int*);
    void Cblacs_gridinit(int*, const char*, int, int);
    void Cblacs_pcoord(int, int, int*, int*);
    void Cblacs_gridexit(int);
    void Cblacs_barrier(int, const char*);
    void Cdgerv2d(int, int, int, double*, int, int, int);
    void Cdgesd2d(int, int, int, double*, int, int, int);
    void Cblacs_gridmap(int*, int*, int, int, int);
    void Cblacs_exit(int);    
    int numroc_(int*, int*, int*, int*, int*);
}

int main(int argc, char **argv)
{

    int     context, desca[9], descb[9], ib, info, ja, laf, lda, ldb,
        lwork, mb, mype, n, nb, npcol, npe, nprow, nrhs;
    double  b[4], d[4], dl[4], du[4];
    double  *af, *work;
    char    trans = 'N';
    /* Set array dimensions and blocking */
    n = 8;            /* dimension of the problem */
    lda = 4;          /* leading dimension of A */
    ldb = 4;          /* leading dimension of B */
    nrhs = 1;         /* number of right-hand sides  */
    npcol = 2;        /* number of processor columns */
    ja = 1;           /* offset for A */
    ib = ja;          /* offset for B */
    mb = 4;           /* blocking */
    nb = 4;           /* blocking */
    laf = 12 * npcol + 3 * nb;
    af = (double *)malloc(laf*sizeof(double));
    lwork = 10 * npcol + 4 * nrhs;
    work = (double *)malloc(lwork*sizeof(double));
    /* Start BLACS */

    Cblacs_pinfo(&mype, &npe);
    Cblacs_get(0, 0, &context);
    Cblacs_gridinit(&context, "R", 1, npe);
    if (mype == 0){
        /* PE = 0 gets D(1:4), DL(1:4), DU(1:4) and B(1:4) */
        d[0] = 1.8180; d[1] = 1.6602; d[2] = 1.3420; d[3] = 1.2897;
        dl[0] = 0.0000; dl[1] = 0.8385; dl[2] = 0.5681; dl[3] = 0.3704;
        du[0] = 0.6946; du[1] = 0.4449; du[2] = 0.5466; du[3] = 0.7027;
        b[0] = 1.0; b[1] = 2.0; b[2] = 3.0; b[3] = 4.0;
    }
    else if (mype == 1){
        /* PE = 1 gets D(5:8), DL(5:8), DU(5:8) and B(5:8) */
        d[0] = 1.3412; d[1] = 1.5341; d[2] = 1.7271; d[3] = 1.3093;
        dl[0] = 0.7027; dl[1] = 0.5466; dl[2] = 0.4449; dl[3] = 0.6946;
        du[0] = 0.3704; du[1] = 0.5681; du[2] = 0.8385; du[3] = 0.0000;
        b[0] = 5.0; b[1] = 6.0; b[2] = 7.0; b[3] = 8.0;
    }
    /* Array descriptor for A (D, DL and DU) */
    desca[0] = 501; desca[1] = context; desca[2] = n; desca[3] = nb;
    desca[4] = 0; desca[5] = lda; desca[6] = 0;
    /* Array descriptor for B */
    descb[0] = 502; descb[1] = context; descb[2] = n; descb[3] = nb;
    descb[4] = 0; descb[5] = ldb; descb[6] = 0;
    /* Factorization */
    pddttrf(&n, dl, d, du, &ja, desca, af, &laf, work, &lwork, &info);
    /* Solution */
    //pddttrs(&trans, &n, &nrhs, dl, d, du, &ja, desca, b, &ib, descb,
    //  af, &laf, work, &lwork, &info);
    printf("MYPE=%i: x[:] = %7.4f %7.4f %7.4f %7.4fn",
        mype, b[0], b[1], b[2], b[3]);
    Cblacs_gridexit(context);
    Cblacs_exit(0);
}

来自英特尔论坛,只是一个评论(评论区太大了),希望它能帮助将来的人。

[yhu5@prc-mic01 scalapack]$run -np 2 ./a.out
MYPE=0: x[:] =  1.0000  2.0000  3.0000  4.0000
MYPE=1: x[:] =  5.0000  6.0000  7.0000  8.0000

可以运行2个秩。按照代码设计。

如果以1 rank运行,它将返回错误。

[yhu5@prc-mic01 scalapack]$ mpirun -np 1 ./a.out
{    0,    0}:  On entry to
PDDTTRF, D&C alg.: only 1 block per proc parameter number    1 had an illegal value
pddttrf problem! Info -1
MYPE=0: x[:] =  1.0000  2.0000  3.0000  4.0000

似乎代码一开始出错了Cblacs_get(0, 0, &context);,因为他是{ -1, -1}