#include #include #include #include "mpi.h" #include "mkl.h" #include "mkl_cluster_sparse_solver.h" // mpiicc -g -DMKL_ILP64 -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -lmkl_blacs_intelmpi_ilp64 -liomp5 -lpthread -lm -ldl cluster_solver_calculix_simple.c // mpirun -check-mpi -np 2 ./a.out void dummy_cluster_sparse_solver(); int main (void) { /* -------------------------------------------------------------------- */ /* .. Init MPI. */ /* -------------------------------------------------------------------- */ /* Auxiliary variables. */ int mpi_stat = 0; int argc = 0; int comm, rank; char** argv; mpi_stat = MPI_Init( &argc, &argv ); mpi_stat = MPI_Comm_rank( MPI_COMM_WORLD, &rank ); comm = MPI_Comm_c2f( MPI_COMM_WORLD ); //if ( rank > 0 ) { dummy_cluster_sparse_solver(); } /* Matrix data. */ MKL_INT n = 8; MKL_INT ia[9] = { 0, 4, 7, 9, 11, 14, 16, 17, 18 }; MKL_INT ja[18] = { 0, 2, 5, 6, /* index of non-zeros in 0 row*/ 1, 2, 4, /* index of non-zeros in 1 row*/ 2, 7, /* index of non-zeros in 2 row*/ 3, 6, /* index of non-zeros in 3 row*/ 4, 5, 6, /* index of non-zeros in 4 row*/ 5, 7, /* index of non-zeros in 5 row*/ 6, /* index of non-zeros in 6 row*/ 7 /* index of non-zeros in 7 row*/ }; float a[18] = { 7.0, /*0*/ 1.0, /*0*/ /*0*/ 2.0, 7.0, /*0*/ -4.0, 8.0, /*0*/ 2.0, /*0*/ /*0*/ /*0*/ 1.0, /*0*/ /*0*/ /*0*/ /*0*/ 5.0, 7.0, /*0*/ /*0*/ 9.0, /*0*/ 5.0, 1.0, 5.0, /*0*/ -1.0, /*0*/ 5.0, 11.0, /*0*/ 5.0 }; MKL_INT mtype = -2; /* set matrix type to "real symmetric indefinite matrix" */ MKL_INT nrhs = 1; /* number of right hand sides. */ float b[8], x[8], bs[8], res, res0; /* RHS and solution vectors. */ /* Internal solver memory pointer pt * 32-bit: int pt[64] or void *pt[64]; * 64-bit: long int pt[64] or void *pt[64]; */ void *pt[64] = { 0 }; void *pt_2[64] = { 0 }; /* Cluster Sparse Solver control parameters. */ MKL_INT iparm[64] = { 0 }; MKL_INT maxfct, mnum, phase, msglvl, error1, error2; /* Auxiliary variables. */ float ddum; /* float dummy */ MKL_INT idum; /* Integer dummy. */ MKL_INT i, j; /* -------------------------------------------------------------------- */ /* .. Setup Cluster Sparse Solver control parameters. */ /* -------------------------------------------------------------------- */ iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */ iparm[ 1] = 2; /* Use METIS for fill-in reordering */ iparm[ 5] = 0; /* Write solution into x */ iparm[ 7] = 2; /* Max number of iterative refinement steps */ iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ iparm[10] = 0; /* Don't use nonsymmetric permutation and scaling MPS */ iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */ iparm[18] = -1; /* Output: Mflops for LU factorization */ iparm[27] = 1; /* Single precision mode of Cluster Sparse Solver */ iparm[26] = 1; //gf iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */ iparm[39] = 0; /* Input: matrix/rhs/solution stored on master */ maxfct = 1; /* Maximum number of numerical factorizations. */ mnum = 1; /* Which factorization to use. */ msglvl = 1; /* Print statistical information in file */ error1 = 0; /* Initialize error flag */ error2 = 0; /* Initialize error flag */ /* -------------------------------------------------------------------- */ /* .. Reordering and Symbolic Factorization. This step also allocates */ /* all memory that is necessary for the factorization. */ /* -------------------------------------------------------------------- */ phase = 12; MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); cluster_sparse_solver ( pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &comm, &error1 ); MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); cluster_sparse_solver ( pt_2, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &comm, &error2 ); if ( error1 != 0 ) { printf ("\nERROR during symbolic factorization, pt: %lli", (long long int)error1); mpi_stat = MPI_Finalize(); return 1; } if ( error2 != 0 ) { printf ("\nERROR during symbolic factorization, pt_2: %lli", (long long int)error2); mpi_stat = MPI_Finalize(); return 1; } printf ("\nReordering completed ... "); /* -------------------------------------------------------------------- */ /* .. Back substitution and iterative refinement. */ /* -------------------------------------------------------------------- */ /* Set right hand side to one. */ for ( i = 0; i < n; i++ ) { b[i] = 1.0; x[i] = 0.0; } printf ("\nSolving system..."); phase = 33; MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); cluster_sparse_solver ( pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &comm, &error1 ); if ( error1 != 0 ) { printf ("\nERROR during solution, pt: %lli", (long long int)error1); mpi_stat = MPI_Finalize(); return 4; } printf ("\nThe solution of the system is: "); for ( j = 0; j < n ; j++ ) { printf ( "\n x1[%lli] = % f", (long long int)j, x[j] ); } /* -------------------------------------------------------------------- */ /* .. Termination and release of memory. */ /* -------------------------------------------------------------------- */ phase = -1; /* Release internal memory. */ MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); cluster_sparse_solver ( pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &comm, &error1 ); if ( error1 != 0 ) { printf ("\nERROR during release memory, pt: %lli", (long long int)error1); goto final; } /* -------------------------------------------------------------------- */ /* .. Repeat phase 33 for second matrix */ /* -------------------------------------------------------------------- */ /* phase = 12; MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); cluster_sparse_solver ( pt_2, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &comm, &error2 ); */ phase=33; MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); cluster_sparse_solver ( pt_2, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &comm, &error2 ); printf ("\nThe solution of the system is: "); for ( j = 0; j < n ; j++ ) { printf ( "\n x2[%lli] = % f", (long long int)j, x[j] ); } phase = 1; /* Release internal memory. */ MPI_Send(&phase, 1, MPI_LONG_LONG, 1, 0 , MPI_COMM_WORLD); /* -------------------------------------------------------------------- */ final: if ( error1 != NULL || error2 != NULL ) { printf("\n TEST FAILED\n"); } else { printf("\n TEST PASSED\n"); } mpi_stat = MPI_Finalize(); return error1; } //////////////////////////////////// void dummy_cluster_sparse_solver() { int mpi_stat = 0; int argc = 0; int comm, rank; char** argv; mpi_stat = MPI_Comm_rank( MPI_COMM_WORLD, &rank ); comm = MPI_Comm_c2f( MPI_COMM_WORLD ); /* Matrix data. */ MKL_INT n; MKL_INT *ia; MKL_INT *ja; MKL_INT mtype; MKL_INT nrhs; double *a, *b, *x; void *pt[64] = { 0 }; /* Cluster Sparse Solver control parameters. */ MKL_INT iparm[64] = { 0 }; MKL_INT maxfct, mnum, msglvl, error; double ddum; /* float dummy */ MKL_INT idum; /* Integer dummy. */ MKL_INT phase; MPI_Recv(&phase, 1, MPI_LONG_LONG, 0 , 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE ); while(( phase != 1 )){ printf ( "\nEntering phase %i while loop\n", phase ); cluster_sparse_solver ( pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum,&ddum, &comm, &error ); MPI_Recv(&phase, 1, MPI_LONG_LONG, 0 , 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE ); } // end while mpi_stat = MPI_Finalize(); exit(0); } // end function