SLEPc eigensolver that uses minimal memory and finds ALL eigenvalues of a real symmetric sparse matrix in reasonable time
Hi all !! The reason I have started using SLEPc/PETSc is that I want to reduce the memory requirements of a program that I already have in Intel MKL. The program's purpose is to calculate ALL eigenvalues and eigenvectors of a real symmetric sparse matrix. Intel MKL doesnt have a sparse solver for this (cant solve for eigenvalues of a matrix represented in sparse format), so I am using SLEPc, as it can do this. (NOTE: by "sparsity" I mean ZEROES*100/(ZEROES+NONZEROES). Also the tolerance I am using is 10^-6) Till now I have been using the LAPACK method (EPSLAPACK) for the slepc eigensolver, as I found out that it was the fastest (I configured PETSc to use Intel MKL BLAS/LAPACK routines). But I just found out that it is using a lot of memory irrespective of the sparsity of the matrix. For example, on a matrix of order 9600x9600 no matter how sparse/dense it is (I tested from 0.02% to 75% sparsity), it used 1.3gb of peak memory every time, whereas the part of the program in which the matrix is assembled in AIJ format (in case of 75% sparsity) takes only 550mb. The MKL program takes 3.1gb. Although memory requirement is less, the difference is constant irrespective of sparsity.So,* it seems that the LAPACK method will use memory dependent on the order of the matrix, irrespective of the sparsity.* So I guess its not very suitable for my purpose ? I also tried changing the NCV and MPD (Maximum projected Dimension) parameters hoping they would somehow reduce memory requirements even if at the cost of more time, but whatever I pass NCV and MPD (in EPSSetDimensions) as, they are not considered at all. They remain the same as the default values (as I verified with EPSGetDimensions). So I tried other solvers like the default, krylovschur, lanczos, but I think they are built to find only a few eigenvalues and not all, because they almost never converged to find ALL eigenvalues, even if I gave them a lot of time (big max number of iterations parameter). But their memory requirements seem to be pretty less. For example, I ran the program with EPSLANCZOS on the same 9600x9600 matrix (75% sparsity) with maximum number of iterations as 500,000. It ran it for 6.5 hrs (then I killed it), and at the end it was using 530mb only (although I dont know how near it was to completion). I am somewhat illiterate about the mathematics behind all this, and how the lapack,krylovschur, etc. work and I have no idea if what I am experiencing was expected/obvious, although I half-expected some things by reading the slepc documentation, where a table was given comparing different eigensolvers. Also, I have observed that as the sparsity of the matrix increases (more zeroes), more eigenvalues also become zeroes. Is there a way to exploit this using probably krylovschur/lanczos so that there is no need to converge to all eigenvalues (meaning that the time they will take will be reasonable), but to only the number of eigenvalues that are non-zero,so that I can assume the non-converged eigenvalues to be simply zeroes (but is there a way to know beforehand by sparsity approximately how many non-zeroe eigenvalues will be there, and the corresponding value of maximum iterations that are required ?) Is there a way to do what I want to do through SLEPc ??? * I want ALL eigenvalues and want that the memory requirements of the eigensolver should be less/in accordance with the sparsity of the matrix*, otherwise the very purpose of using a sparse matrix representation is somewhat defeated. I am doing all this in FORTRAN by the way. Thanks a lot in advance !!! Shitij
El 03/08/2011, a las 08:30, Shitij Bhargava escribió:
Hi all !!
The reason I have started using SLEPc/PETSc is that I want to reduce the memory requirements of a program that I already have in Intel MKL. The program's purpose is to calculate ALL eigenvalues and eigenvectors of a real symmetric sparse matrix. Intel MKL doesnt have a sparse solver for this (cant solve for eigenvalues of a matrix represented in sparse format), so I am using SLEPc, as it can do this. (NOTE: by "sparsity" I mean ZEROES*100/(ZEROES+NONZEROES). Also the tolerance I am using is 10^-6)
If you really want to compute ALL eigenvalues, then SLEPc is not the way to go. SLEPc is intended for computing a few eigenpairs, and in some cases a few percentage, e.g. 20% at most. Computing more than that with SLEPc is probably going to be less efficient than other alternatives.
Till now I have been using the LAPACK method (EPSLAPACK) for the slepc eigensolver, as I found out that it was the fastest (I configured PETSc to use Intel MKL BLAS/LAPACK routines). But I just found out that it is using a lot of memory irrespective of the sparsity of the matrix. For example, on a matrix of order 9600x9600 no matter how sparse/dense it is (I tested from 0.02% to 75% sparsity), it used 1.3gb of peak memory every time, whereas the part of the program in which the matrix is assembled in AIJ format (in case of 75% sparsity) takes only 550mb. The MKL program takes 3.1gb. Although memory requirement is less, the difference is constant irrespective of sparsity.So, it seems that the LAPACK method will use memory dependent on the order of the matrix, irrespective of the sparsity. So I guess its not very suitable for my purpose ? I also tried changing the NCV and MPD (Maximum projected Dimension) parameters hoping they would somehow reduce memory requirements even if at the cost of more time, but whatever I pass NCV and MPD (in EPSSetDimensions) as, they are not considered at all. They remain the same as the default values (as I verified with EPSGetDimensions).
The EPSLAPACK solver converts the matrix to dense format. As explained in the documentation, this solver must be used only on small matrices for debugging purposes. NCV and MPD are not relevant in this setting.
So I tried other solvers like the default, krylovschur, lanczos, but I think they are built to find only a few eigenvalues and not all, because they almost never converged to find ALL eigenvalues, even if I gave them a lot of time (big max number of iterations parameter). But their memory requirements seem to be pretty less. For example, I ran the program with EPSLANCZOS on the same 9600x9600 matrix (75% sparsity) with maximum number of iterations as 500,000. It ran it for 6.5 hrs (then I killed it), and at the end it was using 530mb only (although I dont know how near it was to completion).
Krylov-Schur is almost always preferred over Lanczos, since the latter implements an explicit restart which is much worse than implicit restart in Krylov-Schur. Anyway, in both of them if you set NEV=9600 you are building a projected problem (dense matrix) of order equal to the original matrix so you have the same cost as with EPSLAPACK together with the cost of building the Krylov subspace. Again: do not use SLEPc for all eigenvalues!
I am somewhat illiterate about the mathematics behind all this, and how the lapack,krylovschur, etc. work and I have no idea if what I am experiencing was expected/obvious, although I half-expected some things by reading the slepc documentation, where a table was given comparing different eigensolvers.
A general recommendation is to try to understand a little bit how the algorithms work before attempting to use them.
Also, I have observed that as the sparsity of the matrix increases (more zeroes), more eigenvalues also become zeroes. Is there a way to exploit this using probably krylovschur/lanczos so that there is no need to converge to all eigenvalues (meaning that the time they will take will be reasonable), but to only the number of eigenvalues that are non-zero,so that I can assume the non-converged eigenvalues to be simply zeroes (but is there a way to know beforehand by sparsity approximately how many non-zeroe eigenvalues will be there, and the corresponding value of maximum iterations that are required ?)
There is no direct connection between the number of zero entries and the number of zero eigenvalues, unless you have rows/columns with all zeros, in which case the (permuted) matrix can be written as A=[A11 0; 0 0] so it is enough to compute eigenvalues of A11.
Is there a way to do what I want to do through SLEPc ??? I want ALL eigenvalues and want that the memory requirements of the eigensolver should be less/in accordance with the sparsity of the matrix, otherwise the very purpose of using a sparse matrix representation is somewhat defeated. I am doing all this in FORTRAN by the way.
Ask yourself this question: is it really necessary to compute all eigenvalues? If the answer is yes then maybe ScaLAPACK is what you need. Jose
Thanks a lot in advance !!!
Shitij
Thank you very much for your reply. I am a summer trainee, and I talked to my guide about this, and now what I want might be a little different. I really do want all the eigenvalues. There cannot be any compromise on that. But please note that time is not much of a concern here, the much bigger concern is of memory. Even if I can distribute the memory somehow for the computation of all eigenvalues, that is fine. (doesnt matter if it takes a huge amount of time, he is willing to wait). Can any slepc Eigensolver distribute memory this way while solving for all eigenvalues ? I am sorry if what I am asking is obvious, I dont know much about distributed memory, etc. also. I am stupid at this moment, but I am learning. Also, when I ran the slepc program on a small cluster, and then ran the top command, I saw that while solving for eigenvalues, the program was already using upto 500% CPU. Does this mean that the eigensolver was "distributing memory" automatically among different nodes ? I understand that it must be distributing computational effort, but how do I know if it was distributing memory also ? It simply showed me the RES memory usage to be 1.3 gb, and the %MEM to be about 17%. Also, I did not use any MPIXXX data structures. I simply used SeqAIJ. Will using MPIAIJ "distribute memory" while solving for all eigenvalues ? (I understand that it will distribute memory to store the matrix, but the memory bottleneck is eigenvalue solver, so its more important for memory to be distributed then ) My guide said that one node will not have enough memory, but all nodes combined will, and hence the requirement of memory to be distributed. I read about ScaLAPACK, and I have already asked on their forums to confirm if it is suitable for me. I am waiting for an answer. Thank you once again. Shitij On 3 August 2011 14:31, Jose E. Roman <[email protected]> wrote:
El 03/08/2011, a las 08:30, Shitij Bhargava escribió:
Hi all !!
The reason I have started using SLEPc/PETSc is that I want to reduce the memory requirements of a program that I already have in Intel MKL. The program's purpose is to calculate ALL eigenvalues and eigenvectors of a real symmetric sparse matrix. Intel MKL doesnt have a sparse solver for this (cant solve for eigenvalues of a matrix represented in sparse format), so I am using SLEPc, as it can do this. (NOTE: by "sparsity" I mean ZEROES*100/(ZEROES+NONZEROES). Also the tolerance I am using is 10^-6)
If you really want to compute ALL eigenvalues, then SLEPc is not the way to go. SLEPc is intended for computing a few eigenpairs, and in some cases a few percentage, e.g. 20% at most. Computing more than that with SLEPc is probably going to be less efficient than other alternatives.
Till now I have been using the LAPACK method (EPSLAPACK) for the slepc
eigensolver, as I found out that it was the fastest (I configured PETSc to use Intel MKL BLAS/LAPACK routines). But I just found out that it is using a lot of memory irrespective of the sparsity of the matrix. For example, on a matrix of order 9600x9600 no matter how sparse/dense it is (I tested from 0.02% to 75% sparsity), it used 1.3gb of peak memory every time, whereas the part of the program in which the matrix is assembled in AIJ format (in case of 75% sparsity) takes only 550mb. The MKL program takes 3.1gb. Although memory requirement is less, the difference is constant irrespective of sparsity.So, it seems that the LAPACK method will use memory dependent on the order of the matrix, irrespective of the sparsity. So I guess its not very suitable for my purpose ? I also tried changing the NCV and MPD (Maximum projected Dimension) parameters hoping they would somehow reduce memory requirements even if at the cost of more time, but whatever I pass NCV and MPD (in EPSSetDimensions) as, they are not considered at all. They remain the same as the default values (as I verified with EPSGetDimensions).
The EPSLAPACK solver converts the matrix to dense format. As explained in the documentation, this solver must be used only on small matrices for debugging purposes. NCV and MPD are not relevant in this setting.
So I tried other solvers like the default, krylovschur, lanczos, but I
think they are built to find only a few eigenvalues and not all, because they almost never converged to find ALL eigenvalues, even if I gave them a lot of time (big max number of iterations parameter). But their memory requirements seem to be pretty less. For example, I ran the program with EPSLANCZOS on the same 9600x9600 matrix (75% sparsity) with maximum number of iterations as 500,000. It ran it for 6.5 hrs (then I killed it), and at the end it was using 530mb only (although I dont know how near it was to completion).
Krylov-Schur is almost always preferred over Lanczos, since the latter implements an explicit restart which is much worse than implicit restart in Krylov-Schur. Anyway, in both of them if you set NEV=9600 you are building a projected problem (dense matrix) of order equal to the original matrix so you have the same cost as with EPSLAPACK together with the cost of building the Krylov subspace. Again: do not use SLEPc for all eigenvalues!
I am somewhat illiterate about the mathematics behind all this, and how
the lapack,krylovschur, etc. work and I have no idea if what I am experiencing was expected/obvious, although I half-expected some things by reading the slepc documentation, where a table was given comparing different eigensolvers.
A general recommendation is to try to understand a little bit how the algorithms work before attempting to use them.
Also, I have observed that as the sparsity of the matrix increases (more
zeroes), more eigenvalues also become zeroes. Is there a way to exploit this using probably krylovschur/lanczos so that there is no need to converge to all eigenvalues (meaning that the time they will take will be reasonable), but to only the number of eigenvalues that are non-zero,so that I can assume the non-converged eigenvalues to be simply zeroes (but is there a way to know beforehand by sparsity approximately how many non-zeroe eigenvalues will be there, and the corresponding value of maximum iterations that are required ?)
There is no direct connection between the number of zero entries and the number of zero eigenvalues, unless you have rows/columns with all zeros, in which case the (permuted) matrix can be written as A=[A11 0; 0 0] so it is enough to compute eigenvalues of A11.
Is there a way to do what I want to do through SLEPc ??? I want ALL
eigenvalues and want that the memory requirements of the eigensolver should be less/in accordance with the sparsity of the matrix, otherwise the very purpose of using a sparse matrix representation is somewhat defeated. I am doing all this in FORTRAN by the way.
Ask yourself this question: is it really necessary to compute all eigenvalues? If the answer is yes then maybe ScaLAPACK is what you need.
Jose
Thanks a lot in advance !!!
Shitij
El 04/08/2011, a las 08:48, Shitij Bhargava escribió:
Thank you very much for your reply.
I am a summer trainee, and I talked to my guide about this, and now what I want might be a little different.
I really do want all the eigenvalues. There cannot be any compromise on that. But please note that time is not much of a concern here, the much bigger concern is of memory. Even if I can distribute the memory somehow for the computation of all eigenvalues, that is fine. (doesnt matter if it takes a huge amount of time, he is willing to wait). Can any slepc Eigensolver distribute memory this way while solving for all eigenvalues ?
I am sorry if what I am asking is obvious, I dont know much about distributed memory, etc. also. I am stupid at this moment, but I am learning.
Also, when I ran the slepc program on a small cluster, and then ran the top command, I saw that while solving for eigenvalues, the program was already using upto 500% CPU. Does this mean that the eigensolver was "distributing memory" automatically among different nodes ? I understand that it must be distributing computational effort, but how do I know if it was distributing memory also ? It simply showed me the RES memory usage to be 1.3 gb, and the %MEM to be about 17%. Also, I did not use any MPIXXX data structures. I simply used SeqAIJ. Will using MPIAIJ "distribute memory" while solving for all eigenvalues ? (I understand that it will distribute memory to store the matrix, but the memory bottleneck is eigenvalue solver, so its more important for memory to be distributed then ) My guide said that one node will not have enough memory, but all nodes combined will, and hence the requirement of memory to be distributed.
I read about ScaLAPACK, and I have already asked on their forums to confirm if it is suitable for me. I am waiting for an answer.
Thank you once again.
Shitij
Try something like this: $ mpirun -np 8 ./program -eps_nev 4800 -eps_mpd 400 -eps_largest_real $ mpirun -np 8 ./program -eps_nev 4800 -eps_mpd 400 -eps_smallest_real But this may not work for the smallest ones. Jose
Thank you very much for your reply. Unfortunately, I do not have access to the cluster at this moment because there is some problem with it. But I made a seperate program that reads a matrix of size 2500x2500 from a file and calculates eigenvalues.The matrix itself is sparse, symmetric and has non-zero values randomly distributed. I read more about the krylovschur and other solvers, and parallel computing and mpi, and I probably understand what you mean. The eigensolvers are efficient for calculating eigenvalues at the peripheries, and especially quickly (less number of iterations) on the larger side than the smaller side, right ? Although I still dont understand what Maximum Projected Dimension exactly means. I played with mpd values nonetheless, and it doesnt seem to have much effect on memory usage. I am running it on my laptop. But the main problem seems to be that in the program, different threads dont seem to 'cooperate'. Every thread is calculating its own solution and storing its own matrix, as if running its own version of EPSSolve. This is because in the output to top command, each process (I used 2) uses nearly equal amount of memory, which is equal to the memory usage if I run the program with only 1 process. (so I suppose memory distribution is not happening ?) I am using MPIAIJ now. I am putting the code below: (I understand why the program prints things multiple times, but I know I can fix that by something like "if(processID=0) printSomething"....what I am not able to understand is how to distribute memory while solving for eigenvalues, which is the call to EPSSolve. There doesnt seem to be any parameter to tell it to solve while cooperating with other threads) *#include "slepceps.h" #include <stdio.h> #include <time.h> #include <math.h> #define TOL 0.000001 #define MAXIT 1000 int main(int argc,char**argv) { FILE *fp; double entry; int N; int nzval=0,Istart,Iend; Mat A; PetscErrorCode ierr; double *value,kr,ki; int i,j,nconv,nev,maxit,its; int col_nz;/*number of non-zero columns in a row*/ Vec xr, xi; double tol,error,re,im; EPS eps; /*Eigen solver context*/ const EPSType type; int ncv,mpd; /*PetscInt rows;*/ int* columns; SlepcInitialize(&argc,&argv,(char*)0,PETSC_NULL); /*Read matrix from file*/ fp=fopen("sparseMatrix.txt","r"); fscanf(fp,"%d",&N); fscanf(fp,"%lf",&entry);/*not required, but only to move file pointer forward*/ /*Allocate memory to arrays*/ value=(double*)malloc(N*sizeof(double)); columns=(int*)malloc(N*sizeof(int)); ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); ierr=MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); //ierr=MatSeqAIJSetPreallocation(A,N,PETSC_NULL);CHKERRQ(ierr); ierr=MatMPIAIJSetPreallocation(A,N,PETSC_NULL,N,PETSC_NULL);CHKERRQ(ierr);//?????!!!!!!!!!!!!!!! ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); for(i=0;i<N;i++) { col_nz=0; if(i<Iend && i>=Istart) //important, see where interval closed and open { for(j=0;j<N;j++) { fscanf(fp,"%lf",&entry); if(entry!=0) { value[col_nz]=entry; columns[col_nz]=j; col_nz++; } } ierr = MatSetValues(A,1,&i,col_nz,columns,value,INSERT_VALUES);CHKERRQ(ierr); } } fclose(fp); MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); printf("MATRIX ASSMEBLY DONE !!!!!!!! \n\n"); /*Solve the Eigen Problem*/ /***********START TIME************/ clock_t start=clock(); /********************************/ ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr); ierr = EPSSetOperators(eps,A,PETSC_NULL);CHKERRQ(ierr); ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); ierr=EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);CHKERRQ(ierr); ierr=EPSSetDimensions(eps,0.5*N,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);/*N specifies that find ALL eigenvalues*/ ierr = EPSSetTolerances(eps,TOL,500);CHKERRQ(ierr); //ierr = MatDestroy(A);CHKERRQ(ierr); //dont need this anymore ierr = EPSSolve(eps);CHKERRQ(ierr); printf("Solved for half the largest values\n\n\n\n"); ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);CHKERRQ(ierr); ierr = EPSDestroy(eps);CHKERRQ(ierr); ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr); ierr = EPSSetOperators(eps,A,PETSC_NULL);CHKERRQ(ierr); ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); ierr=EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr); ierr=EPSSetDimensions(eps,0.5*N,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);/*N specifies that find ALL eigenvalues*/ ierr = EPSSetTolerances(eps,TOL,500);CHKERRQ(ierr); ierr = EPSSolve(eps);CHKERRQ(ierr); ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);CHKERRQ(ierr); ierr = EPSDestroy(eps);CHKERRQ(ierr); /**********STOP TIME**************/ clock_t end=clock(); /********************************/ ierr = PetscPrintf(PETSC_COMM_WORLD,"Time taken to compute all eigenvalues of a %d x %d matrix is : %lf seconds\n\n",N,N,(double)(end-start)/CLOCKS_PER_SEC );CHKERRQ(ierr); ierr = SlepcFinalize();CHKERRQ(ierr); return 0; } *I* *ran it with: mpirun -np 2 ./slepcEigenMPI -eps_monitor I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension. This is the output I get: (part of the output) *MATRIX ASSMEBLY DONE !!!!!!!! MATRIX ASSMEBLY DONE !!!!!!!! 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 3 EPS nconv=386 first unconverged value (error) 2.81217e-25 (6.07438011e-04) 3 EPS nconv=386 first unconverged value (error) 2.81217e-25 (6.07438011e-04) 4 EPS nconv=420 first unconverged value (error) 2.75201e-27 (1.65345126e-03) 4 EPS nconv=420 first unconverged value (error) 2.75201e-27 (1.65345126e-03) 6 EPS nconv=1001 first unconverged value (error) 1.93852e-27 (8.40346506e-02) 6 EPS nconv=1001 first unconverged value (error) 1.93852e-27 (8.40346506e-02) 7 EPS nconv=1027 first unconverged value (error) 2.00333e-27 (7.40100945e-03) 7 EPS nconv=1027 first unconverged value (error) 2.00333e-27 (7.40100945e-03)* How do I make "different threads store the matrix and more importantly, run EPSSolve together" ? Thats what I want right ? Thank you very much ! Shitij ** On 4 August 2011 16:57, Jose E. Roman <[email protected]> wrote:
El 04/08/2011, a las 08:48, Shitij Bhargava escribió:
Thank you very much for your reply.
I am a summer trainee, and I talked to my guide about this, and now what I want might be a little different.
I really do want all the eigenvalues. There cannot be any compromise on that. But please note that time is not much of a concern here, the much bigger concern is of memory. Even if I can distribute the memory somehow for the computation of all eigenvalues, that is fine. (doesnt matter if it takes a huge amount of time, he is willing to wait). Can any slepc Eigensolver distribute memory this way while solving for all eigenvalues ?
I am sorry if what I am asking is obvious, I dont know much about distributed memory, etc. also. I am stupid at this moment, but I am learning.
Also, when I ran the slepc program on a small cluster, and then ran the top command, I saw that while solving for eigenvalues, the program was already using upto 500% CPU. Does this mean that the eigensolver was "distributing memory" automatically among different nodes ? I understand that it must be distributing computational effort, but how do I know if it was distributing memory also ? It simply showed me the RES memory usage to be 1.3 gb, and the %MEM to be about 17%. Also, I did not use any MPIXXX data structures. I simply used SeqAIJ. Will using MPIAIJ "distribute memory" while solving for all eigenvalues ? (I understand that it will distribute memory to store the matrix, but the memory bottleneck is eigenvalue solver, so its more important for memory to be distributed then ) My guide said that one node will not have enough memory, but all nodes combined will, and hence the requirement of memory to be distributed.
I read about ScaLAPACK, and I have already asked on their forums to confirm if it is suitable for me. I am waiting for an answer.
Thank you once again.
Shitij
Try something like this:
$ mpirun -np 8 ./program -eps_nev 4800 -eps_mpd 400 -eps_largest_real $ mpirun -np 8 ./program -eps_nev 4800 -eps_mpd 400 -eps_smallest_real
But this may not work for the smallest ones.
Jose
On Mon, Aug 8, 2011 at 00:29, Shitij Bhargava <[email protected]> wrote:
I* *ran it with:
mpirun -np 2 ./slepcEigenMPI -eps_monitor
I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension.
This is the output I get: (part of the output) *MATRIX ASSMEBLY DONE !!!!!!!!
MATRIX ASSMEBLY DONE !!!!!!!!
1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04)*
The most likely case is that you have more than one MPI implementation installed and that you are running with a different implementation than you built with. Compare the outputs: $ ldd ./slepcEigenMPI $ which mpirun
Thank you Jed. That was indeed the problem. I installed a separate MPI for PETSc/SLEPc, but was running my program with a default, already installed one. Now, I have a different question. What I want to do is this: 1. Only 1 process, say root, calculates the matrix in SeqAIJ format 2. Then root creates the EPS context, eps and initializes,sets parameters, problem type,etc. properly 3. After this the root process broadcasts this eps object to other processes 4. I use EPSSolve to solve for eigenvalues (all process together in cooperation resulting in memory distribution) 5. I get the results from root is this possible ? I am not able to broadcast the EPS object, because it is not an MPI_DataType. Is there any PETSc/SLEPc function for this ? I am avoiding using MPIAIJ because that will mean making many changes in the existing code, including the numerous write(*,*) statements (i would have to convert them to PetscPrint in FORTRAN or something like that). So I want a single process to handle matrix generation and assembly, but want to solve the eigenproblem in parallel by different processes. Running the subroutine EPSSolve in parallel and hence distribute memory is the only reason why I want to use MPI. Thanks a lot !! Shitij On 8 August 2011 11:05, Jed Brown <[email protected]> wrote:
On Mon, Aug 8, 2011 at 00:29, Shitij Bhargava <[email protected]>wrote:
I* *ran it with:
mpirun -np 2 ./slepcEigenMPI -eps_monitor
I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension.
This is the output I get: (part of the output) *MATRIX ASSMEBLY DONE !!!!!!!!
MATRIX ASSMEBLY DONE !!!!!!!!
1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04)*
The most likely case is that you have more than one MPI implementation installed and that you are running with a different implementation than you built with. Compare the outputs:
$ ldd ./slepcEigenMPI $ which mpirun
El 08/08/2011, a las 09:14, Shitij Bhargava escribió:
Thank you Jed. That was indeed the problem. I installed a separate MPI for PETSc/SLEPc, but was running my program with a default, already installed one.
Now, I have a different question. What I want to do is this:
1. Only 1 process, say root, calculates the matrix in SeqAIJ format 2. Then root creates the EPS context, eps and initializes,sets parameters, problem type,etc. properly 3. After this the root process broadcasts this eps object to other processes 4. I use EPSSolve to solve for eigenvalues (all process together in cooperation resulting in memory distribution) 5. I get the results from root
is this possible ? I am not able to broadcast the EPS object, because it is not an MPI_DataType. Is there any PETSc/SLEPc function for this ? I am avoiding using MPIAIJ because that will mean making many changes in the existing code, including the numerous write(*,*) statements (i would have to convert them to PetscPrint in FORTRAN or something like that). So I want a single process to handle matrix generation and assembly, but want to solve the eigenproblem in parallel by different processes. Running the subroutine EPSSolve in parallel and hence distribute memory is the only reason why I want to use MPI.
Thanks a lot !!
Shitij
No, you cannot use a solver in parallel with a matrix in SeqAIJ format. The matrix must be MPIAIJ. If you want to generate the matrix only in process 0, you can do MPI_Comm_rank(PETSC_COMM_WORLD,&rank) and then enclose the matrix generation (setvalues only) in an if clause: if (!rank) { ierr = MatSetValues(A,...);CHKERRQ(ierr); } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); However, it is better to do the matrix generation in parallel also. Jose
On Aug 8, 2011, at 2:14 AM, Shitij Bhargava wrote:
Thank you Jed. That was indeed the problem. I installed a separate MPI for PETSc/SLEPc, but was running my program with a default, already installed one.
Now, I have a different question. What I want to do is this:
1. Only 1 process, say root, calculates the matrix in SeqAIJ format 2. Then root creates the EPS context, eps and initializes,sets parameters, problem type,etc. properly 3. After this the root process broadcasts this eps object to other processes 4. I use EPSSolve to solve for eigenvalues (all process together in cooperation resulting in memory distribution) 5. I get the results from root
We do have an undocumented routine MatDistribute_MPIAIJ(MPI_Comm comm,Mat gmat,PetscInt m,MatReuse reuse,Mat *inmat) in src/mat/impls/aij/mpi/mpiaij.c that will take a SeqAIJ matrix and distribute it over a larger MPI communicator. Note that you cannot create the EPS context etc on a the root process and then broadcast the object but once the matrix is distributed you can simple create the EPS context etc on the parallel communicator where the matrix is and run with that. Barry
is this possible ? I am not able to broadcast the EPS object, because it is not an MPI_DataType. Is there any PETSc/SLEPc function for this ? I am avoiding using MPIAIJ because that will mean making many changes in the existing code, including the numerous write(*,*) statements (i would have to convert them to PetscPrint in FORTRAN or something like that). So I want a single process to handle matrix generation and assembly, but want to solve the eigenproblem in parallel by different processes. Running the subroutine EPSSolve in parallel and hence distribute memory is the only reason why I want to use MPI.
Thanks a lot !!
Shitij
On 8 August 2011 11:05, Jed Brown <[email protected]> wrote: On Mon, Aug 8, 2011 at 00:29, Shitij Bhargava <[email protected]> wrote: I ran it with:
mpirun -np 2 ./slepcEigenMPI -eps_monitor
I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension.
This is the output I get: (part of the output) MATRIX ASSMEBLY DONE !!!!!!!!
MATRIX ASSMEBLY DONE !!!!!!!!
1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04)
The most likely case is that you have more than one MPI implementation installed and that you are running with a different implementation than you built with. Compare the outputs:
$ ldd ./slepcEigenMPI $ which mpirun
Thanks Jose, Barry. I tried what you said, but that gives me an error: *[0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: Can only get local values, trying 9!* This is probably because here I am trying to insert all rows of the matrix through process 0, but process 0 doesnt own all the rows. In any case, this seems very "unnatural", so I am using MPIAIJ the right way as you said, where I assemble the MPIAIJ matrix in parallel instead of only on one process. I have done that actually, and am running the code on the cluster right now. Its going to take a long long time to finish, so I cant confirm some of my doubts, which I am asking below: 1. If I run the code with 1 process, and say it takes M memory (peak) while solving for eigenvalues, then when I run it with N processes, each will take nearly M/N memory (peak) (probably a little more) right ? And for doing this, I dont have to use any special MPI stuff....the fact that I am using MPIAIJ, and building the EPS object from it, and then calling EPSSolve() is enough ? I mean EPSSolve() is internally in some way distributing memory and computation effort automatically when I use MPIAIJ, and run the code with many processes, right ? This confusion is there because when I use top, while running the code with 8 processes, each of them showed me nearly 250 mb initially, but each has grown to use 270 mb in about 70 minutes. I understand that the method krylovschur is such that memory requirements increase slowly, but the peak on any process will be less (than if I ran only one process), right ? (Even though their memory requirements are growing, they will grow to some M/N only, right ?) Actually the fact that in this case, each of the process creates its own EPS context, initializes it itself, and then calls EPSSolve() itself without any "interaction" with other processes makes me wonder if they really are working together, or just individually (I would have verified this myself, but the program will take way too much time, and I know I would have to kill it sooner or later).....or the fact that they initialize their own EPS context with THEIR part of the MPI is enough to make them "cooperate and work together" ? (Although I think this is what Barry meant in that last post, but I am not too sure) I am not too comfortable with the MPI way of thinking right now, probably this is why I have this confusion. Anyways, I cant thank you guys enough. I would have been scrounging through documentation again and again to no avail if you guys had not helped me the way you did. The responses were always prompt, always to the point (even though my questions were sometimes not, probably because I didnt completely understand the problems I was facing.....but you always knew what I was asking) and very clear. At this moment, I dont know much about PETSc/SLEPc myself, but I will be sure to contribute back to this list when I do. I have nothing but sincere gratitude for you guys. Thank you very much ! Shitij On 9 August 2011 00:58, Barry Smith <[email protected]> wrote:
On Aug 8, 2011, at 2:14 AM, Shitij Bhargava wrote:
Thank you Jed. That was indeed the problem. I installed a separate MPI for PETSc/SLEPc, but was running my program with a default, already installed one.
Now, I have a different question. What I want to do is this:
1. Only 1 process, say root, calculates the matrix in SeqAIJ format 2. Then root creates the EPS context, eps and initializes,sets parameters, problem type,etc. properly 3. After this the root process broadcasts this eps object to other processes 4. I use EPSSolve to solve for eigenvalues (all process together in cooperation resulting in memory distribution) 5. I get the results from root
We do have an undocumented routine MatDistribute_MPIAIJ(MPI_Comm comm,Mat gmat,PetscInt m,MatReuse reuse,Mat *inmat) in src/mat/impls/aij/mpi/mpiaij.c that will take a SeqAIJ matrix and distribute it over a larger MPI communicator.
Note that you cannot create the EPS context etc on a the root process and then broadcast the object but once the matrix is distributed you can simple create the EPS context etc on the parallel communicator where the matrix is and run with that.
Barry
is this possible ? I am not able to broadcast the EPS object, because it
is not an MPI_DataType. Is there any PETSc/SLEPc function for this ? I am avoiding using MPIAIJ because that will mean making many changes in the existing code, including the numerous write(*,*) statements (i would have to convert them to PetscPrint in FORTRAN or something like that).
So I want a single process to handle matrix generation and assembly, but want to solve the eigenproblem in parallel by different processes. Running the subroutine EPSSolve in parallel and hence distribute memory is the only reason why I want to use MPI.
Thanks a lot !!
Shitij
On 8 August 2011 11:05, Jed Brown <[email protected]> wrote: On Mon, Aug 8, 2011 at 00:29, Shitij Bhargava <[email protected]> wrote: I ran it with:
mpirun -np 2 ./slepcEigenMPI -eps_monitor
I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension.
This is the output I get: (part of the output) MATRIX ASSMEBLY DONE !!!!!!!!
MATRIX ASSMEBLY DONE !!!!!!!!
1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04)
The most likely case is that you have more than one MPI implementation installed and that you are running with a different implementation than you built with. Compare the outputs:
$ ldd ./slepcEigenMPI $ which mpirun
El 09/08/2011, a las 09:54, Shitij Bhargava escribió:
Thanks Jose, Barry.
I tried what you said, but that gives me an error:
[0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: Can only get local values, trying 9!
This is probably because here I am trying to insert all rows of the matrix through process 0, but process 0 doesnt own all the rows.
In any case, this seems very "unnatural", so I am using MPIAIJ the right way as you said, where I assemble the MPIAIJ matrix in parallel instead of only on one process. I have done that actually, and am running the code on the cluster right now. Its going to take a long long time to finish, so I cant confirm some of my doubts, which I am asking below:
1. If I run the code with 1 process, and say it takes M memory (peak) while solving for eigenvalues, then when I run it with N processes, each will take nearly M/N memory (peak) (probably a little more) right ? And for doing this, I dont have to use any special MPI stuff....the fact that I am using MPIAIJ, and building the EPS object from it, and then calling EPSSolve() is enough ? I mean EPSSolve() is internally in some way distributing memory and computation effort automatically when I use MPIAIJ, and run the code with many processes, right ? This confusion is there because when I use top, while running the code with 8 processes, each of them showed me nearly 250 mb initially, but each has grown to use 270 mb in about 70 minutes. I understand that the method krylovschur is such that memory requirements increase slowly, but the peak on any process will be less (than if I ran only one process), right ? (Even though their memory requirements are growing, they will grow to some M/N only, right ?)
The solver allocates some dynamic memory when the actual computation starts, so it is normal that you see a growth in the memory footprint. No further increase should be observed afterwards. Jose
Actually the fact that in this case, each of the process creates its own EPS context, initializes it itself, and then calls EPSSolve() itself without any "interaction" with other processes makes me wonder if they really are working together, or just individually (I would have verified this myself, but the program will take way too much time, and I know I would have to kill it sooner or later).....or the fact that they initialize their own EPS context with THEIR part of the MPI is enough to make them "cooperate and work together" ? (Although I think this is what Barry meant in that last post, but I am not too sure)
I am not too comfortable with the MPI way of thinking right now, probably this is why I have this confusion.
Anyways, I cant thank you guys enough. I would have been scrounging through documentation again and again to no avail if you guys had not helped me the way you did. The responses were always prompt, always to the point (even though my questions were sometimes not, probably because I didnt completely understand the problems I was facing.....but you always knew what I was asking) and very clear. At this moment, I dont know much about PETSc/SLEPc myself, but I will be sure to contribute back to this list when I do. I have nothing but sincere gratitude for you guys.
Thank you very much !
Shitij
On Aug 9, 2011, at 2:54 AM, Shitij Bhargava wrote:
Thanks Jose, Barry.
I tried what you said, but that gives me an error:
[0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: Can only get local values, trying 9!
This is probably because here I am trying to insert all rows of the matrix through process 0, but process 0 doesnt own all the rows.
In any case, this seems very "unnatural", so I am using MPIAIJ the right way as you said, where I assemble the MPIAIJ matrix in parallel instead of only on one process. I have done that actually, and am running the code on the cluster right now. Its going to take a long long time to finish,
It shouldn't take a long time to finish. Are you sure you are creating all the objects with the PETSC_COMM_WORLD and not PETSC_COMM_SELF? Have you done the correct matrix preallocation http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#efficient-assem... Is each process generating just its part of the matrix? Barry
so I cant confirm some of my doubts, which I am asking below:
1. If I run the code with 1 process, and say it takes M memory (peak) while solving for eigenvalues, then when I run it with N processes, each will take nearly M/N memory (peak) (probably a little more) right ? And for doing this, I dont have to use any special MPI stuff....the fact that I am using MPIAIJ, and building the EPS object from it, and then calling EPSSolve() is enough ? I mean EPSSolve() is internally in some way distributing memory and computation effort automatically when I use MPIAIJ, and run the code with many processes, right ? This confusion is there because when I use top, while running the code with 8 processes, each of them showed me nearly 250 mb initially, but each has grown to use 270 mb in about 70 minutes. I understand that the method krylovschur is such that memory requirements increase slowly, but the peak on any process will be less (than if I ran only one process), right ? (Even though their memory requirements are growing, they will grow to some M/N only, right ?)
Actually the fact that in this case, each of the process creates its own EPS context, initializes it itself, and then calls EPSSolve() itself without any "interaction" with other processes makes me wonder if they really are working together, or just individually (I would have verified this myself, but the program will take way too much time, and I know I would have to kill it sooner or later).....or the fact that they initialize their own EPS context with THEIR part of the MPI is enough to make them "cooperate and work together" ? (Although I think this is what Barry meant in that last post, but I am not too sure)
I am not too comfortable with the MPI way of thinking right now, probably this is why I have this confusion.
Anyways, I cant thank you guys enough. I would have been scrounging through documentation again and again to no avail if you guys had not helped me the way you did. The responses were always prompt, always to the point (even though my questions were sometimes not, probably because I didnt completely understand the problems I was facing.....but you always knew what I was asking) and very clear. At this moment, I dont know much about PETSc/SLEPc myself, but I will be sure to contribute back to this list when I do. I have nothing but sincere gratitude for you guys.
Thank you very much !
Shitij
On 9 August 2011 00:58, Barry Smith <[email protected]> wrote:
On Aug 8, 2011, at 2:14 AM, Shitij Bhargava wrote:
Thank you Jed. That was indeed the problem. I installed a separate MPI for PETSc/SLEPc, but was running my program with a default, already installed one.
Now, I have a different question. What I want to do is this:
1. Only 1 process, say root, calculates the matrix in SeqAIJ format 2. Then root creates the EPS context, eps and initializes,sets parameters, problem type,etc. properly 3. After this the root process broadcasts this eps object to other processes 4. I use EPSSolve to solve for eigenvalues (all process together in cooperation resulting in memory distribution) 5. I get the results from root
We do have an undocumented routine MatDistribute_MPIAIJ(MPI_Comm comm,Mat gmat,PetscInt m,MatReuse reuse,Mat *inmat) in src/mat/impls/aij/mpi/mpiaij.c that will take a SeqAIJ matrix and distribute it over a larger MPI communicator.
Note that you cannot create the EPS context etc on a the root process and then broadcast the object but once the matrix is distributed you can simple create the EPS context etc on the parallel communicator where the matrix is and run with that.
Barry
is this possible ? I am not able to broadcast the EPS object, because it is not an MPI_DataType. Is there any PETSc/SLEPc function for this ? I am avoiding using MPIAIJ because that will mean making many changes in the existing code, including the numerous write(*,*) statements (i would have to convert them to PetscPrint in FORTRAN or something like that). So I want a single process to handle matrix generation and assembly, but want to solve the eigenproblem in parallel by different processes. Running the subroutine EPSSolve in parallel and hence distribute memory is the only reason why I want to use MPI.
Thanks a lot !!
Shitij
On 8 August 2011 11:05, Jed Brown <[email protected]> wrote: On Mon, Aug 8, 2011 at 00:29, Shitij Bhargava <[email protected]> wrote: I ran it with:
mpirun -np 2 ./slepcEigenMPI -eps_monitor
I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension.
This is the output I get: (part of the output) MATRIX ASSMEBLY DONE !!!!!!!!
MATRIX ASSMEBLY DONE !!!!!!!!
1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04)
The most likely case is that you have more than one MPI implementation installed and that you are running with a different implementation than you built with. Compare the outputs:
$ ldd ./slepcEigenMPI $ which mpirun
Thank you both for your replies. Actually, Barry, you are right. On eight CPUs, as it turns out, it will not take as long as I imagined. For a 9600x9600 matrix, to solve for half the largest eigenvalues, it took about 300 minutes. Although it would have taken much more time than this for solving half the smallest eigenvalues (I had to kill it at total time of 600 minutes). This is still much much longer than the LAPACK method, which takes (for calculating all the eigenvalues at once) about 90 minutes (at the cost of much memory, which cant even be distributed -- which is unacceptable. Also, I suppose you probably didnt know that I have to calculate ALL the eigenvalues). But still, I was expecting that it wouldnt complete even in "days", that is why I said it would take a very long time....I forgot to take into account the fact that it was now running nearly 8 times faster. (in one earlier instance I had to kill the process when I was running the same program on one process only (was not using MPI), as it had run over 13 hours and still hadnt completed !!) And yes, I am making all the objects with PETSC_COMM_WORLD only, also, the top command shows me eight processes running at nearly 100% CPU during call to EPSSolve(). (I suppose that verifies what I am saying ? Though I am not sure...) I have done preallocation carefully, and verified that that part doesnt take any extra time...so, actually when I referred to "time" here, I was talking about time taken by the eigensolver only (because I was already satisfied by the time taken for matrix generation and assembly) Each process is not "generating" its part of the matrix only at this moment. I am doing things in a very lousy way right now. Each process is generating all the rows of the matrix, but inserts only the rows which fall into its ownership range....but as I said, matrix generation and assembly is not the time bottleneck in the program....eigensolver is....I am satisfied with the time matrix generation and assembly is taking, although I understand that ideally matrix generation should also be done in parallel. My time for this summer training has ended, so I'll probably fix that later. The amount of work that I have been able to accomplish in this training, has been good because of your help. I might just have been stuck in some error otherwise. Thank you very much once again !! Shitij On 10 August 2011 02:02, Barry Smith <[email protected]> wrote:
On Aug 9, 2011, at 2:54 AM, Shitij Bhargava wrote:
Thanks Jose, Barry.
I tried what you said, but that gives me an error:
[0]PETSC ERROR: --------------------- Error Message
[0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: Can only get local values, trying 9!
This is probably because here I am trying to insert all rows of the matrix through process 0, but process 0 doesnt own all the rows.
In any case, this seems very "unnatural", so I am using MPIAIJ the right way as you said, where I assemble the MPIAIJ matrix in parallel instead of only on one process. I have done that actually, and am running the code on the cluster right now. Its going to take a long long time to finish,
It shouldn't take a long time to finish. Are you sure you are creating all the objects with the PETSC_COMM_WORLD and not PETSC_COMM_SELF? Have you done the correct matrix preallocation http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#efficient-assem... Is each process generating just its part of the matrix?
Barry
so I cant confirm some of my doubts, which I am asking below:
1. If I run the code with 1 process, and say it takes M memory (peak) while solving for eigenvalues, then when I run it with N processes, each will take nearly M/N memory (peak) (probably a little more) right ? And for doing this, I dont have to use any special MPI stuff....the fact that I am using MPIAIJ, and building the EPS object from it, and then calling EPSSolve() is enough ? I mean EPSSolve() is internally in some way distributing memory and computation effort automatically when I use MPIAIJ, and run the code with many processes, right ? This confusion is there because when I use top, while running the code with 8 processes, each of them showed me nearly 250 mb initially, but each has grown to use 270 mb in about 70 minutes. I understand that the method krylovschur is such that memory requirements increase slowly, but the peak on any process will be less (than if I ran only one process), right ? (Even though their memory requirements are growing, they will grow to some M/N only, right ?)
Actually the fact that in this case, each of the process creates its own EPS context, initializes it itself, and then calls EPSSolve() itself without any "interaction" with other processes makes me wonder if they really are working together, or just individually (I would have verified this myself, but the program will take way too much time, and I know I would have to kill it sooner or later).....or the fact that they initialize their own EPS context with THEIR part of the MPI is enough to make them "cooperate and work together" ? (Although I think this is what Barry meant in that last post, but I am not too sure)
I am not too comfortable with the MPI way of thinking right now, probably this is why I have this confusion.
Anyways, I cant thank you guys enough. I would have been scrounging through documentation again and again to no avail if you guys had not helped me the way you did. The responses were always prompt, always to the point (even though my questions were sometimes not, probably because I didnt completely understand the problems I was facing.....but you always knew what I was asking) and very clear. At this moment, I dont know much about PETSc/SLEPc myself, but I will be sure to contribute back to this list when I do. I have nothing but sincere gratitude for you guys.
Thank you very much !
Shitij
On 9 August 2011 00:58, Barry Smith <[email protected]> wrote:
On Aug 8, 2011, at 2:14 AM, Shitij Bhargava wrote:
Thank you Jed. That was indeed the problem. I installed a separate MPI for PETSc/SLEPc, but was running my program with a default, already installed one.
Now, I have a different question. What I want to do is this:
1. Only 1 process, say root, calculates the matrix in SeqAIJ format 2. Then root creates the EPS context, eps and initializes,sets parameters, problem type,etc. properly 3. After this the root process broadcasts this eps object to other processes 4. I use EPSSolve to solve for eigenvalues (all process together in cooperation resulting in memory distribution) 5. I get the results from root
We do have an undocumented routine MatDistribute_MPIAIJ(MPI_Comm comm,Mat gmat,PetscInt m,MatReuse reuse,Mat *inmat) in src/mat/impls/aij/mpi/mpiaij.c that will take a SeqAIJ matrix and distribute it over a larger MPI communicator.
Note that you cannot create the EPS context etc on a the root process and then broadcast the object but once the matrix is distributed you can simple create the EPS context etc on the parallel communicator where the matrix is and run with that.
Barry
is this possible ? I am not able to broadcast the EPS object, because
it is not an MPI_DataType. Is there any PETSc/SLEPc function for this ? I am avoiding using MPIAIJ because that will mean making many changes in the existing code, including the numerous write(*,*) statements (i would have to convert them to PetscPrint in FORTRAN or something like that).
So I want a single process to handle matrix generation and assembly, but want to solve the eigenproblem in parallel by different processes. Running the subroutine EPSSolve in parallel and hence distribute memory is the only reason why I want to use MPI.
Thanks a lot !!
Shitij
On 8 August 2011 11:05, Jed Brown <[email protected]> wrote: On Mon, Aug 8, 2011 at 00:29, Shitij Bhargava <[email protected]> wrote: I ran it with:
mpirun -np 2 ./slepcEigenMPI -eps_monitor
I didnt do exactly what you said, because the matrix generation part in the actual program is quite time consuming itself. But I assume what I am doing is equivalent to what you meant to do? Also, I put MPD as PETSC_DECIDE, because I didnt know what to put it for this matrix dimension.
This is the output I get: (part of the output) MATRIX ASSMEBLY DONE !!!!!!!!
MATRIX ASSMEBLY DONE !!!!!!!!
1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 1 EPS nconv=98 first unconverged value (error) 1490.88 (1.73958730e-05) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04) 2 EPS nconv=282 first unconverged value (error) 3.04636e-27 (2.49532175e-04)
The most likely case is that you have more than one MPI implementation installed and that you are running with a different implementation than you built with. Compare the outputs:
$ ldd ./slepcEigenMPI $ which mpirun
On Thu, Aug 11, 2011 at 02:23, Shitij Bhargava <[email protected]> wrote:
On eight CPUs, as it turns out, it will not take as long as I imagined. For a 9600x9600 matrix, to solve for half the largest eigenvalues, it took about 300 minutes. Although it would have taken much more time than this for solving half the smallest eigenvalues (I had to kill it at total time of 600 minutes). This is still much much longer than the LAPACK method, which takes (for calculating all the eigenvalues at once) about 90 minutes (at the cost of much memory, which cant even be distributed -- which is unacceptable.
You should try Elemental for distributed-memory dense eigensolvers. http://code.google.com/p/elemental/
Also, I suppose you probably didnt know that I have to calculate ALL the eigenvalues).
Why do you need all of them? What underlying problem are you solving?
participants (4)
-
Barry Smith -
Jed Brown -
Jose E. Roman -
Shitij Bhargava