5
Я пишу небольшой тестовый код для диагонализации параллельной матрицы с использованием алгоритма разделения и покоя ScaLAPACK PDSYEVD на C. Однако я новичок в ScaLAPACK и, глядя на источник, кажется довольно страшным количеством переменные, для которых я не мог найти хороших примеров. Тип матриц, которые мне нужно диагонализировать, являются реальными, симметричными, довольно разреженными и размером ~ 1000х1000.Диагональ матрицы скайпака (pdsyevd)
Простой (последовательный) 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? По-видимому, это частичная матрица собственного вектора, но почему это не написано в A, как в DSYEV?
- У всех процессов должен быть одинаковый размер блока? Может ли быть перекрытие подматриц или что я делаю для матриц простого размера, например?
- Что находится в 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;
}
ИМО этот вопрос будет лучше всего разделить на отдельные вопросы для отдельных частей. –