我正在使用 C 中的 ScaLAPACK 的分而治之算法PDSYEVD编写一个用于并行矩阵对角化的小型测试代码。但是我是 ScaLAPACK 的新手，并且查看源代码似乎需要设置的变量数量相当可怕，但我无法找到好的任何例子。我需要对角化的矩阵类型是真实的、对称的、相当稀疏的，并且大小约为 1000x1000。
一个简单的（串行）基于 LAPACK 的代码如下所示：
/* "matrix diagonalization using lapack" */
#include <stdio.h>
#include <stdlib.h>
/* DSYEV prototype */
extern void dsyev_(char* jobz, char* uplo, int* n, double* a, int* lda,
double* w, double* work, int* lwork, int* info );
#define n 4
int main(int argc, char *argv[])
{
int i,j,info,lda,lwork,N;
double A[n][n], a[n*n], work[100*n],w[n];
/* initialize matrices */
for(i=0;i<n;i++) { for(j=0;j<n;j++) A[i][j] = (i+j); }
/* diagonalize */
N=n; lda=n; lwork=10*n; info=0;
for(i=0;i<n;i++) { for(j=0;j<n;j++) a[i+j*n]=A[i][j]; }
dsyev_("V","L",&N,a,&lda,w,work,&lwork,&info);
/* print result */
for(i=0;i<n;i++) {
printf("w[%d] = %8.4f | v[%d] = [ ",i,w[i],i);
for(j=0;j<n;j++) printf("%6.2f ",a[j+i*n]);
printf("]\n");
}
return 0;
}
从那里我想去 pdsyevd，它应该被称为：
SUBROUTINE PDSYEVD( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO )
很多事情我都不清楚。对这些问题中的任何一个或几个问题的答案将不胜感激。
- 有没有人碰巧有一个例子？
- 在 ScaLAPACK 文档中它说：“[...] PDSYEVD 假设同质系统 [...] 异构系统可能会返回不正确的结果而没有任何错误消息。”。在这种情况下，“同质/异质”系统是什么意思？
- 在文档的末尾它说：“分布式子矩阵必须验证一些对齐属性，即 [...]”。这个表达式中的术语甚至不在输入中，我怎么知道这是否满足？MB_A/MB_Z？
- 如何为每个流程选择子矩阵？有什么道理吗？我想它们的大小应该或多或少相等，但是行/列/正方形/...？对于更稀疏的区域可能更大......？
- Z/IZ/IJ/DESCZ 在输入中做了什么？我已经为每个进程提供了一个子矩阵 A/IA/JA/DESCA 来处理，那么为什么是这个 Z？它似乎是部分特征向量矩阵，但为什么不像 DSYEV 中那样写入 A？
- 所有进程都必须具有相同的块大小吗？子矩阵是否可以重叠，或者我对素数大小的矩阵该怎么做？
- DESCA中有什么？另外，什么是DLEN？它似乎是 A 的输入“数组描述符”。我不知道该怎么做。
到目前为止，或多或少是我所拥有的，并没有做太多：
/* "matrix diagonalization using scalapack" */
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* PDSYEVD prototype */
extern void pdsyevd_( char* jobz, char* uplo, int* n, double* a, int* ia, int* ja, int* desca,
double* w, double* z, int* iz, int* jz, int* descz, double* work,
int* lwork, double* iwork, int* liwork, int* info );
#define n 4
int main(int argc, char *argv[])
{
int numprocs,myid,i,j,k,index,info,lwork,liwork,N;
int ia,ja,desca[n];
int iz,jz,descz[n];
double A[n][n],w[n];
double *a,*z,*work,*iwork;
/* set up MPI stuff */
MPI_Init(&argc,&argv); // initialize MPI
MPI_Comm_size(MPI_COMM_WORLD,&numprocs); // find out size of world
MPI_Comm_rank(MPI_COMM_WORLD,&myid); // determine current mpi proc id
/* initialize matrices */
for(i=0;i<n;i++) { for(j=0;j<n;j++) A[i][j] = (i+j); }
/* determine submatrices */
/* let's use square blocks, of equal size:
there number of processors (mp^2) should be
- square (mp^2 = 1,2,4,9,16,..)
- with mp being a divisor of n
In case we have 4 procs and a 4x4 matrix for example
(A11 A12 A13 A14)
(A21 A22 A23 A24)
A = (A31 A32 A33 A34)
(A41 A42 A43 A44)
( a1 a2 )
= ( a3 a3 )
where (A11 A12)
a1 = (A21 A22)
is the submatrix sent to process 1 for example
in this case :
n=4
nprocs=4 --> mp=2
len = 2
*/
int mp = sqrt(numprocs);
if(numprocs!=mp*mp) {printf("ERROR: numprocs (%d) should be square.\n",numprocs); return 1; }
if(n%mp!=0) {printf("ERROR: mp (%d) should be divisor of n (%d).\n",mp,n); return 1; }
int len = n/mp;
a = malloc((n*n+1)*sizeof(double));
z = malloc((n*n+1)*sizeof(double));
//a = malloc((len*len+1)*sizeof(double));
//z = malloc((len*len+1)*sizeof(double));
ia=(myid*len)%n; // from 0 to n in steps of len
ja=(int)(myid/mp)*len;
printf("proc %d has a submatrix of size %d, with starting indices (%d,%d)\n",myid,len,ia,ja);
iz=ia;
jz=ja;
for(i=ia;i<ia+len;i++) { for(j=ja,k=0;j<ja+len;j++,k++) a[k]=A[i][j]; }
for(i=iz;i<iz+len;i++) { for(j=jz,k=0;j<jz+len;j++,k++) z[k]=A[i][j]; }
/* diagonalize */
N=n; lwork=n*n; liwork=n*n; info=0;
work = malloc(lwork*sizeof(double));
iwork = malloc(liwork*sizeof(double));
for(i=0;i<n;i++) { for(j=0;j<n;j++) a[i+j*n]=A[i][j]; }
pdsyevd_("V","U",&len,a,&ia,&ja,desca,w,z,&iz,&jz,descz,work,&lwork,iwork,&liwork,&info);
if(!myid)printf("info = %d\n",info);
/* gather & print results */
double wf[n];
MPI_Reduce(w,wf,n,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
if(myid==0) {
for(i=0;i<n;i++) {
printf("wf[%d] = %8.4f | v[%d] = [ ",i,wf[i],i);
for(j=0;j<n;j++) printf("%6.2f ",a[j+i*n]);
printf("]\n");
}
}
free(a); free(z); free(work); free(iwork);
/* clean up */
MPI_Finalize(); /* MPI Programs end with MPI Finalize; this is a weak synchronization point */
return 0;
}