1 /* Modified by Fernando Lorenzo on December 13, 2003 to compile with Superlu 3.0 */ 2 /* ++++++ This file will not work if Superlu 1.0 is used ++++++++++++*/ 3 4 /* 5 Copyright (C) 2000 Dennis Roddeman 6 email: dennis.roddeman@feat.nl 7 8 This program is free software; you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation; either version 2 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program; if not, write to the Free Software Foundation 20 59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA 21 */ 22 23 #include "tnsuplu.h" 24 25 #if ( SUPERLU_USE || SUPERLU_MT_USE || SUPERLU_DIST_USE ) 26 27 #if SUPERLU_USE 28 #include "slu_ddefs.h" 29 #elif SUPERLU_MT_USE 30 #include "pdsp_defs.h" 31 #elif SUPERLU_DIST_USE 32 #include <math.h> 33 #include "superlu_ddefs.h" 34 #endif 35 #endif 36 #include "slu_util.h" 37 38 long int solve_superlu( double *superlu_A, int *superlu_asub, 39 int *superlu_xa, int superlu_nnz, double solve_b[], long int solve_nlocal, 40 long int nthread ) 41 42 { 43 int n ,info; 44 45 #if SUPERLU_DIST_USE 46 47 SuperMatrix A; 48 int_t nprow=1, npcol=1, nrhs=1; 49 SuperLUStat_t stat; 50 superlu_options_t options; 51 ScalePermstruct_t ScalePermstruct; 52 LUstruct_t LUstruct; 53 gridinfo_t grid; 54 double *berr; 55 int iam, count; 56 #else 57 58 SuperMatrix A, L, U, B; 59 NCformat *Astore; 60 superlu_options_t options; 61 SuperLUStat_t stat; 62 mem_usage_t mem_usage; 63 int *perm_r, *perm_c; 64 int perm_spec; 65 int nrhs; 66 67 #endif 68 69 /* Create Matrix A in the format expected by SuperLU */ 70 #if SUPERLU_DIST_USE 71 72 /* get number of processors nprow, npcol */ 73 MPI_Comm_size(MPI_COMM_WORLD, &count); 74 nprow=count; npcol=1; 75 /* Initialize superlu process grid */ 76 superlu_gridinit(MPI_COMM_WORLD,nprow,npcol,&grid); 77 iam = grid.iam; 78 /* if we dont belong to this grid , bug out */ 79 if (iam >= nprow*npcol ) goto out; 80 81 #endif 82 n = (int) solve_nlocal; 83 #if SUPERLU_DIST_USE 84 85 if ( !iam ) { 86 /* broadcast matrices to other processes */ 87 MPI_Bcast(&n,1,mpi_int_t,0,grid.comm); 88 MPI_Bcast(&superlu_nnz,1,mpi_int_t,0,grid.comm); 89 MPI_Bcast(superlu_A,superlu_nnz,MPI_DOUBLE,0,grid.comm); 90 MPI_Bcast(superlu_asub,superlu_nnz,mpi_int_t,0,grid.comm); 91 MPI_Bcast(superlu_xa,n+1,mpi_int_t,0,grid.comm); 92 } else { 93 /* receive matrices from proc 0 */ 94 MPI_Bcast(&n,1,mpi_int_t,0,grid.comm); 95 MPI_Bcast(&superlu_nnz,1,mpi_int_t,0,grid.comm); 96 /* allocate storage for compressed column rep */ 97 dallocateA(n,superlu_nnz,&superlu_A,&superlu_asub,&superlu_xa); 98 MPI_Bcast(superlu_A,superlu_nnz,MPI_DOUBLE,0,grid.comm); 99 MPI_Bcast(superlu_asub,superlu_nnz,mpi_int_t,0,grid.comm); 100 MPI_Bcast(superlu_xa,n+1,mpi_int_t,0,grid.comm); 101 } 102 103 #endif 104 /* Set the default input options: 105 options.Fact = DOFACT; 106 options.Equil = YES; 107 options.ColPerm = COLAMD; 108 options.DiagPivotThresh = 1.0; 109 options.Trans = NOTRANS; 110 options.IterRefine = NOREFINE; 111 options.SymmetricMode = NO; 112 options.PivotGrowth = NO; 113 options.ConditionNumber = NO; 114 options.PrintStat = YES; 115 */ 116 set_default_options(&options); 117 /* create compressed column matrix for A */ 118 dCreate_CompCol_Matrix( &A, n, n, superlu_nnz, superlu_A, superlu_asub, 119 superlu_xa, SLU_NC, SLU_D, SLU_GE ); 120 /* Astore = A.Store; 121 printf("Dimension %dx%d; # nonzeros %d\n", A.nrow, A.ncol ); */ 122 123 /* Create right-hand side matrix B. */ 124 nrhs = 1; 125 #if ( SUPERLU_USE || SUPERLU_MT_USE ) 126 127 dCreate_Dense_Matrix( &B, n, nrhs, solve_b, n, SLU_DN, SLU_D, SLU_GE ); 128 129 if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[]."); 130 if ( !(perm_r = intMalloc(n)) ) ABORT("Malloc fails for perm_r[]."); 131 132 #endif 133 134 /* ---------- Solve ----------*/ 135 #if SUPERLU_USE 136 137 /* Initialize Stats vars */ 138 StatInit(&stat); 139 140 dgssv( &options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info); 141 142 #elif SUPERLU_MT_USE 143 pdgssv( nthread, &A, perm_c, perm_r, &L, &U, &B, &info ); 144 #elif SUPERLU_DIST_USE 145 /* set the default input options */ 146 set_default_options(&options); 147 /* Initialize ScalePermstruct and LUstruct */ 148 ScalePermstructInit(n,n,&ScalePermstruct); 149 LUstructInit(n,n,&LUstruct); 150 /* Initialize Stats vars */ 151 PStatInit(&stat); 152 /* call the linear eqn solver */ 153 if (!(berr = doubleMalloc(nrhs)) ) ABORT("Malloc fails for berr[]."); 154 pdgssvx_ABglobal(&options, &A, &ScalePermstruct, solve_b, n, nrhs, 155 &grid, &LUstruct, berr, &stat, &info); 156 /* print the stats */ 157 PStatPrint(&stat,&grid); 158 159 #endif 160 161 /* De-allocate storage */ 162 #if ( SUPERLU_USE || SUPERLU_MT_USE ) 163 164 SUPERLU_FREE(perm_r); 165 SUPERLU_FREE(perm_c); 166 #elif SUPERLU_DIST_USE 167 PStatFree(&stat); 168 Destroy_LU(n,&grid, &LUstruct); 169 ScalePermstructFree(&ScalePermstruct); 170 LUstructFree(&LUstruct); 171 SUPERLU_FREE(berr); 172 173 #endif 174 Destroy_CompCol_Matrix(&A); 175 176 #if ( SUPERLU_USE || SUPERLU_MT_USE ) 177 178 Destroy_SuperMatrix_Store(&B); 179 Destroy_SuperNode_Matrix(&L); 180 Destroy_CompCol_Matrix(&U); 181 182 #endif 183 184 #if SUPERLU_DIST_USE 185 186 /* release the process grid */ 187 out: 188 superlu_gridexit(&grid); 189 190 #endif 191 192 if ( info == 0 ) { 193 194 printf(" Solution Found - L\\U MB %.3f\ttotal MB needed %.3f\n", 195 mem_usage.for_lu/1e6, mem_usage.total_needed/1e6); 196 } else { 197 printf("dgssv() error returns INFO= %d\n", info); 198 } 199 if ( info <= n ) { /* factorization completes */ 200 dQuerySpace(&L, &U, &mem_usage); 201 printf("L\\U MB %.3f\ttotal MB needed %.3f\n", 202 mem_usage.for_lu/1e6, mem_usage.total_needed/1e6); 203 } 204 205 return 1; 206 } 207 208