1C 2C This file is part of MUMPS 5.1.2, released 3C on Mon Oct 2 07:37:01 UTC 2017 4C 5C 6C Copyright 1991-2017 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria, 7C University of Bordeaux. 8C 9C This version of MUMPS is provided to you free of charge. It is 10C released under the CeCILL-C license: 11C http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html 12C 13 SUBROUTINE SMUMPS_SIMSCALEABS(IRN_loc, JCN_loc, A_loc, NZ_loc, 14 & M, N, NUMPROCS, MYID, COMM, 15 & RPARTVEC, CPARTVEC, 16 & RSNDRCVSZ, CSNDRCVSZ, REGISTRE, 17 & IWRK, IWRKSZ, 18 & INTSZ, RESZ, OP, 19 & ROWSCA, COLSCA, WRKRC, ISZWRKRC, 20 & SYM, NB1, NB2, NB3, EPS, 21 & ONENORMERR,INFNORMERR) 22C---------------------------------------------------------------------- 23C IF SYM=0 CALLs unsymmetric variant SMUMPS_SIMSCALEABSUNS. 24C IF SYM=2 CALLS symmetric variant where only one of a_ij and a_ji 25C is stored. SMUMPS_SIMSCALEABSSYM 26C--------------------------------------------------------------------- 27C For details, see the two subroutines below 28C SMUMPS_SIMSCALEABSUNS and SMUMPS_SIMSCALEABSSYM 29C --------------------------------------------------------------------- 30C 31 IMPLICIT NONE 32 INCLUDE 'mpif.h' 33 INTEGER(8) NZ_loc 34 INTEGER IWRKSZ, ISZWRKRC 35 INTEGER M, N, OP 36 INTEGER NUMPROCS, MYID, COMM 37 INTEGER INTSZ, RESZ 38 INTEGER IRN_loc(NZ_loc) 39 INTEGER JCN_loc(NZ_loc) 40 REAL A_loc(NZ_loc) 41 INTEGER RPARTVEC(M) 42 INTEGER RSNDRCVSZ(2*NUMPROCS) 43 INTEGER CPARTVEC(N) 44 INTEGER CSNDRCVSZ(2*NUMPROCS) 45 INTEGER IWRK(IWRKSZ) 46 INTEGER REGISTRE(12) 47 REAL ROWSCA(M) 48 REAL COLSCA(N) 49 REAL WRKRC(ISZWRKRC) 50 REAL ONENORMERR,INFNORMERR 51C LOCALS 52C IMPORTANT POINTERS 53C FOR the scaling phase 54 INTEGER SYM, NB1, NB2, NB3 55 REAL EPS 56C EXTERNALS 57 EXTERNAL SMUMPS_SIMSCALEABSUNS,SMUMPS_SIMSCALEABSSYM, 58 & SMUMPS_INITREAL 59C MUST HAVE IT 60 INTEGER I 61 IF(SYM.EQ.0) THEN 62 CALL SMUMPS_SIMSCALEABSUNS(IRN_loc, JCN_loc, A_loc, 63 & NZ_loc, 64 & M, N, NUMPROCS, MYID, COMM, 65 & RPARTVEC, CPARTVEC, 66 & RSNDRCVSZ, CSNDRCVSZ, REGISTRE, 67 & IWRK, IWRKSZ, 68 & INTSZ, RESZ, OP, 69 & ROWSCA, COLSCA, WRKRC, ISZWRKRC, 70 & NB1, NB2, NB3, EPS, 71 & ONENORMERR, INFNORMERR) 72 ELSE 73 CALL SMUMPS_SIMSCALEABSSYM(IRN_loc, JCN_loc, A_loc, 74 & NZ_loc, 75 & N, NUMPROCS, MYID, COMM, 76 & RPARTVEC, 77 & RSNDRCVSZ, REGISTRE, 78 & IWRK, IWRKSZ, 79 & INTSZ, RESZ, OP, 80 & ROWSCA, WRKRC, ISZWRKRC, 81 & NB1, NB2, NB3, EPS, 82 & ONENORMERR, INFNORMERR) 83 DO I=1,N 84 COLSCA(I) = ROWSCA(I) 85 ENDDO 86 ENDIF 87 RETURN 88 END SUBROUTINE SMUMPS_SIMSCALEABS 89 SUBROUTINE SMUMPS_SIMSCALEABSUNS(IRN_loc, JCN_loc, A_loc, NZ_loc, 90 & M, N, NUMPROCS, MYID, COMM, 91 & RPARTVEC, CPARTVEC, 92 & RSNDRCVSZ, CSNDRCVSZ, REGISTRE, 93 & IWRK, IWRKSZ, 94 & INTSZ, RESZ, OP, 95 & ROWSCA, COLSCA, WRKRC, ISZWRKRC, 96 & NB1, NB2, NB3, EPS, 97 & ONENORMERR, INFNORMERR) 98C---------------------------------------------------------------------- 99C Input parameters: 100C M, N: size of matrix (in general M=N, but the algorithm 101C works for rectangular matrices as well (norms other than 102C inf-norm are not possible mathematically in this case). 103C NUMPROCS, MYID, COMM: guess what are those 104C RPARTVEC: row partvec to be filled when OP=1 105C CPARTVEC: col partvec to be filled when OP=1 106C RSNDRCVSZ: send recv sizes for row operations. 107C to be filled when OP=1 108C CSNDRCVSZ: send recv sizes for col operations. 109C to be filled when OP=1 110C REGISTRE: to store some pointers (size etc) 111C IWRK: working space. when OP=1 IWRKSZ.GE.4*MAXMN 112C when OP=2 INTSZ portion is used. Thus, IWRKSZ>INTSZ 113C when OP=2 114C IWRKSZ: size 115C INTSZ: to be computed when OP=1, necessary integer space to run 116C scaling algo when OP=2 117C RESZ: to be computed when OP=1, necessary real space to run 118C scaling algo when OP=2 119C OP: 120C =1 estimation of memory and construction of partvecs 121C writes into RPARTVEC,CPARTVEC,RSNDRCVSZ,CSNDRCVSZ,REGISTRE 122C does not access WRKRC, uses IWRK as workspace 123C computes INTSZ and RESZ. 124C =2 Compute scalings 125C restores pointers from REGISTRE, 126C stores communication structure in IWRK (from the start). 127C 128C ROWSCA: space for row scaling factor; has size M 129C COLSCA: space for col scaling factor; has size N 130C WRKRC: real working space. when OP=1, is not accessed. Thus, it 131C can be declared to be of size 1 at OP=1 call. 132C ISZWRKRC: size 133C SYM: is matrix symmetric 134C NB1, NB2, NB3: algo runs 135C NB1 iters of inf-norm (default 1/1), 136C NB2 iters of 1-norm (default 3/10), 137C NB3 iters of inf-norm (default 3/10). 138C in succession. 139C EPS: tolerance for concergence. 140C IF EPS < 0.R0 then does not test convergence. 141C If convergence occured during the first set of inf-norm 142C iterations, we start performing one-norm iterations. 143C If convergence occured during the one-norm iterations, 144C we start performing the second set of inf-norm iterations. 145C If convergence occured during the second set of inf-norm, 146C we prepare to return. 147C ONENORMERR : error in one norm scaling (associated with the scaling 148C arrays of the previous iterations), 149C INFNORMERR : error in inf norm scaling (associated with the scaling 150C arrays of the previous iterations). 151C--------------------------------------------------------------------- 152C On input: 153C OP=1==>Requirements 154C IWRKSZ.GE.4*MAXMN 155C RPARTVEC of size M 156C CPARTVEC of size N 157C RSNDRCVSZ of size 2*NUMPROCS 158C CSNDRCVSZ of size 2*NUMPROCS 159C REGISTRE of size 12 160C 161C OP=2==>Requirements 162C INTSZ .GE. REGISTRE(11) 163C RESZ .GE. REGISTRE(12) 164C--------------------------------------------------------------------- 165C On output: 166C ROWSCA and COLSCA 167C at processor 0 of COMM: complete factors. 168C at other processors : only the ROWSCA(i) or COLSCA(j) 169C for which there is a nonzero a_i* or a_*j are useful. 170C ONENORMERR : error in one norm scaling 171C = -1.0 if iter2=0. 172C INFNORMERR : error in inf norm scaling 173C = inf norm error at iter3 if iter3 > 0 174C = inf norm error at iter1 if iter1 > 0, iter3=0 175C = -1.0 if iter1=iter3=0 176C --------------------------------------------------------------------- 177C References: 178C The scaling algorithms are based on those discussed in 179C [1] D. Ruiz, "A scaling algorithm to equilibrate both rows and 180C columns norms in matrices", Tech. Rep. Rutherford 181C Appleton Laboratory, Oxon, UK and ENSEEIHT-IRIT, 182C Toulouse, France, RAL-TR-2001-034 and RT/APO/01/4, 2001. 183C [2] D. Ruiz and B. Ucar, "A symmetry preserving algorithm for 184C matrix scaling", in preparation as of Jan'08. 185C 186C The parallelization approach is discussed in 187C [3] P. R. Amestoy, I. S. Duff, D. Ruiz, and B. Ucar, 188C "A parallel matrix scaling algorithm". 189C In proceedings of VECPAR'08-International Meeting-High 190C Performance Computing for Computational Science, Jan'08. 191C and was supported by ANR-SOLSTICE project (ANR-06-CIS6-010) 192C --------------------------------------------------------------------- 193 IMPLICIT NONE 194 INCLUDE 'mpif.h' 195 INTEGER(8) NZ_loc 196 INTEGER IWRKSZ, INTSZ, ISZWRKRC 197 INTEGER M, N, OP 198 INTEGER NUMPROCS, MYID, COMM 199 INTEGER RESZ 200 INTEGER IRN_loc(NZ_loc) 201 INTEGER JCN_loc(NZ_loc) 202 REAL A_loc(NZ_loc) 203 INTEGER RPARTVEC(M) 204 INTEGER CPARTVEC(N) 205 INTEGER RSNDRCVSZ(2*NUMPROCS) 206 INTEGER CSNDRCVSZ(2*NUMPROCS) 207 INTEGER REGISTRE(12) 208 INTEGER IWRK(IWRKSZ) 209 REAL ROWSCA(M) 210 REAL COLSCA(N) 211 REAL WRKRC(ISZWRKRC) 212 REAL ONENORMERR,INFNORMERR 213C LOCALS 214 INTEGER IRSNDRCVNUM, ORSNDRCVNUM 215 INTEGER ICSNDRCVNUM, OCSNDRCVNUM 216 INTEGER IRSNDRCVVOL, ORSNDRCVVOL 217 INTEGER ICSNDRCVVOL, OCSNDRCVVOL 218 INTEGER INUMMYR, INUMMYC 219C IMPORTANT POINTERS 220 INTEGER IMYRPTR,IMYCPTR 221 INTEGER IRNGHBPRCS, IRSNDRCVIA,IRSNDRCVJA 222 INTEGER ORNGHBPRCS, ORSNDRCVIA,ORSNDRCVJA 223 INTEGER ICNGHBPRCS, ICSNDRCVIA,ICSNDRCVJA 224 INTEGER OCNGHBPRCS, OCSNDRCVIA,OCSNDRCVJA 225 INTEGER ISTATUS, REQUESTS, TMPWORK 226 INTEGER ITDRPTR, ITDCPTR, ISRRPTR 227 INTEGER OSRRPTR, ISRCPTR, OSRCPTR 228C FOR the scaling phase 229 INTEGER NB1, NB2, NB3 230 REAL EPS 231C Iteration vars 232 INTEGER ITER, IR, IC 233 INTEGER(8) :: NZIND 234 REAL ELM 235C COMM TAGS.... 236 INTEGER TAG_COMM_COL 237 PARAMETER(TAG_COMM_COL=100) 238 INTEGER TAG_COMM_ROW 239 PARAMETER(TAG_COMM_ROW=101) 240 INTEGER TAG_ITERS 241 PARAMETER(TAG_ITERS=102) 242C FUNCTIONS 243 EXTERNAL SMUMPS_CREATEPARTVEC, 244 & SMUMPS_NUMVOLSNDRCV, 245 & SMUMPS_SETUPCOMMS, 246 & SMUMPS_FINDNUMMYROWCOL, 247 & SMUMPS_CHKCONVGLO, 248 & SMUMPS_CHK1CONV, 249 & SMUMPS_FILLMYROWCOLINDICES, 250 & SMUMPS_INITREAL, 251 & SMUMPS_INITREALLST, 252 & SMUMPS_DOCOMMINF, 253 & SMUMPS_DOCOMM1N 254 INTEGER SMUMPS_CHKCONVGLO 255 INTEGER SMUMPS_CHK1CONV 256 REAL SMUMPS_ERRSCALOC 257 REAL SMUMPS_ERRSCA1 258 INTRINSIC abs 259 REAL RONE, RZERO 260 PARAMETER(RONE=1.0E0,RZERO=0.0E0) 261C TMP VARS 262 INTEGER RESZR, RESZC 263 INTEGER INTSZR, INTSZC 264 INTEGER MAXMN 265 INTEGER I, IERROR 266 REAL ONEERRROW, ONEERRCOL, ONEERRL, ONEERRG 267 REAL INFERRROW, INFERRCOL, INFERRL, INFERRG 268 INTEGER OORANGEIND 269 INFERRG = -RONE 270 ONEERRG = -RONE 271 OORANGEIND = 0 272 MAXMN = M 273 IF(MAXMN < N) MAXMN = N 274C Create row partvec and col partvec 275 IF(OP == 1) THEN 276 IF(NUMPROCS > 1) THEN 277C Check done outside 278C IF(IWRKSZ.LT.4*MAXMN) THEN ERROR.... 279 CALL SMUMPS_CREATEPARTVEC(MYID, NUMPROCS, COMM, 280 & IRN_loc, JCN_loc, NZ_loc, 281 & RPARTVEC, M, N, 282 & IWRK, IWRKSZ) 283 CALL SMUMPS_CREATEPARTVEC(MYID, NUMPROCS, COMM, 284 & JCN_loc, IRN_loc, NZ_loc, 285 & CPARTVEC, N, M, 286 & IWRK, IWRKSZ) 287C Compute sndrcv sizes, store them for later use 288 CALL SMUMPS_NUMVOLSNDRCV(MYID, NUMPROCS, M, RPARTVEC, 289 & NZ_loc, IRN_loc, N, JCN_loc, 290 & IRSNDRCVNUM,IRSNDRCVVOL, 291 & ORSNDRCVNUM,ORSNDRCVVOL, 292 & IWRK,IWRKSZ, 293 & RSNDRCVSZ(1), RSNDRCVSZ(1+NUMPROCS), COMM) 294 CALL SMUMPS_NUMVOLSNDRCV(MYID, NUMPROCS, N, CPARTVEC, 295 & NZ_loc, JCN_loc, M, IRN_loc, 296 & ICSNDRCVNUM,ICSNDRCVVOL, 297 & OCSNDRCVNUM,OCSNDRCVVOL, 298 & IWRK,IWRKSZ, 299 & CSNDRCVSZ(1), CSNDRCVSZ(1+NUMPROCS), COMM) 300 CALL SMUMPS_FINDNUMMYROWCOL(MYID, NUMPROCS, COMM, 301 & IRN_loc, JCN_loc, NZ_loc, 302 & RPARTVEC, CPARTVEC, M, N, 303 & INUMMYR, 304 & INUMMYC, 305 & IWRK, IWRKSZ) 306 INTSZR = IRSNDRCVNUM + ORSNDRCVNUM + 307 & IRSNDRCVVOL + ORSNDRCVVOL + 308 & 2*(NUMPROCS+1) + INUMMYR 309 INTSZC = ICSNDRCVNUM + OCSNDRCVNUM + 310 & ICSNDRCVVOL + OCSNDRCVVOL + 311 & 2*(NUMPROCS+1) + INUMMYC 312 INTSZ = INTSZR + INTSZC + MAXMN + 313 & (MPI_STATUS_SIZE +1) * NUMPROCS 314 ELSE 315C NUMPROCS IS 1 316 IRSNDRCVNUM = 0 317 ORSNDRCVNUM = 0 318 IRSNDRCVVOL = 0 319 ORSNDRCVVOL = 0 320 INUMMYR = 0 321 ICSNDRCVNUM = 0 322 OCSNDRCVNUM = 0 323 ICSNDRCVVOL = 0 324 OCSNDRCVVOL = 0 325 INUMMYC = 0 326 INTSZ = 0 327 ENDIF 328C CALCULATE NECESSARY REAL SPACE 329 RESZR = M + IRSNDRCVVOL + ORSNDRCVVOL 330 RESZC = N + ICSNDRCVVOL + OCSNDRCVVOL 331 RESZ = RESZR + RESZC 332C CALCULATE NECESSARY INT SPACE 333C The last maxmn is tmpwork for setup comm and fillmyrowcol 334 REGISTRE(1) = IRSNDRCVNUM 335 REGISTRE(2) = ORSNDRCVNUM 336 REGISTRE(3) = IRSNDRCVVOL 337 REGISTRE(4) = ORSNDRCVVOL 338 REGISTRE(5) = ICSNDRCVNUM 339 REGISTRE(6) = OCSNDRCVNUM 340 REGISTRE(7) = ICSNDRCVVOL 341 REGISTRE(8) = OCSNDRCVVOL 342 REGISTRE(9) = INUMMYR 343 REGISTRE(10) = INUMMYC 344 REGISTRE(11) = INTSZ 345 REGISTRE(12) = RESZ 346 ELSE 347C else of op=1. That is op=2 now. 348C restore the numbers 349 IRSNDRCVNUM = REGISTRE(1) 350 ORSNDRCVNUM = REGISTRE(2) 351 IRSNDRCVVOL = REGISTRE(3) 352 ORSNDRCVVOL = REGISTRE(4) 353 ICSNDRCVNUM = REGISTRE(5) 354 OCSNDRCVNUM = REGISTRE(6) 355 ICSNDRCVVOL = REGISTRE(7) 356 OCSNDRCVVOL = REGISTRE(8) 357 INUMMYR = REGISTRE(9) 358 INUMMYC = REGISTRE(10) 359 IF(NUMPROCS > 1) THEN 360C Check done outsize 361C IF(INTSZ < REGISTRE(11)) THEN ERROR 362C IF(RESZ < REGISTRE(12)) THEN ERROR 363C Fill up myrows and my colsX 364 CALL SMUMPS_FILLMYROWCOLINDICES(MYID, NUMPROCS,COMM, 365 & IRN_loc, JCN_loc, NZ_loc, 366 & RPARTVEC, CPARTVEC, M, N, 367 & IWRK(1), INUMMYR, 368 & IWRK(1+INUMMYR), INUMMYC, 369 & IWRK(1+INUMMYR+INUMMYC), IWRKSZ-INUMMYR-INUMMYC ) 370 IMYRPTR = 1 371 IMYCPTR = IMYRPTR + INUMMYR 372C Set up comm and run. 373C set pointers in iwrk (4 parts) 374C 375C ROWS [---------------------------------------------] 376 IRNGHBPRCS = IMYCPTR+ INUMMYC 377 IRSNDRCVIA = IRNGHBPRCS+IRSNDRCVNUM 378 IRSNDRCVJA = IRSNDRCVIA + NUMPROCS+1 379 ORNGHBPRCS = IRSNDRCVJA + IRSNDRCVVOL 380 ORSNDRCVIA = ORNGHBPRCS + ORSNDRCVNUM 381 ORSNDRCVJA = ORSNDRCVIA + NUMPROCS + 1 382C COLS [---------------------------------------------] 383 ICNGHBPRCS = ORSNDRCVJA + ORSNDRCVVOL 384 ICSNDRCVIA = ICNGHBPRCS + ICSNDRCVNUM 385 ICSNDRCVJA = ICSNDRCVIA + NUMPROCS+1 386 OCNGHBPRCS = ICSNDRCVJA + ICSNDRCVVOL 387 OCSNDRCVIA = OCNGHBPRCS + OCSNDRCVNUM 388 OCSNDRCVJA = OCSNDRCVIA + NUMPROCS + 1 389C 390C MPI [-----------------] 391 REQUESTS = OCSNDRCVJA + OCSNDRCVVOL 392 ISTATUS = REQUESTS + NUMPROCS 393C Check done outside 394C IF(ISTATUS + NUMPROCS * MPI_STATUS_SIZE - 1>INTSZ) THEN 395C write(6,*) "Bora: ", ISTATUS + 396C & NUMPROCS * MPI_STATUS_SIZE - 1,INTSZ 397C write(6,*) "Bora : TODO. scimscaent_33 REPORT ERROR" 398C CALL flush(6) 399C ENDIF 400C 401C TMPWRK [-----------------] 402 TMPWORK = ISTATUS + MPI_STATUS_SIZE * NUMPROCS 403 CALL SMUMPS_SETUPCOMMS(MYID, NUMPROCS, M, RPARTVEC, 404 & NZ_loc, IRN_loc,N, JCN_loc, 405 & IRSNDRCVNUM, IRSNDRCVVOL, 406 & IWRK(IRNGHBPRCS),IWRK(IRSNDRCVIA),IWRK(IRSNDRCVJA), 407 & ORSNDRCVNUM, ORSNDRCVVOL, 408 & IWRK(ORNGHBPRCS),IWRK(ORSNDRCVIA),IWRK(ORSNDRCVJA), 409 & RSNDRCVSZ(1), RSNDRCVSZ(1+NUMPROCS), 410 & IWRK(TMPWORK), 411 & IWRK(ISTATUS), IWRK(REQUESTS), 412 & TAG_COMM_ROW, COMM) 413 CALL SMUMPS_SETUPCOMMS(MYID, NUMPROCS, N, CPARTVEC, 414 & NZ_loc, JCN_loc, M, IRN_loc, 415 & ICSNDRCVNUM, ICSNDRCVVOL, 416 & IWRK(ICNGHBPRCS), 417 & IWRK(ICSNDRCVIA), 418 & IWRK(ICSNDRCVJA), 419 & OCSNDRCVNUM, OCSNDRCVVOL, 420 & IWRK(OCNGHBPRCS),IWRK(OCSNDRCVIA),IWRK(OCSNDRCVJA), 421 & CSNDRCVSZ(1), CSNDRCVSZ(1+NUMPROCS), 422 & IWRK(TMPWORK), 423 & IWRK(ISTATUS), IWRK(REQUESTS), 424 & TAG_COMM_COL, COMM) 425 CALL SMUMPS_INITREAL(ROWSCA, M, RZERO) 426 CALL SMUMPS_INITREAL(COLSCA, N, RZERO) 427 CALL SMUMPS_INITREALLST(ROWSCA, M, 428 & IWRK(IMYRPTR),INUMMYR, RONE) 429 CALL SMUMPS_INITREALLST(COLSCA, N, 430 & IWRK(IMYCPTR),INUMMYC, RONE) 431 ELSE 432 CALL SMUMPS_INITREAL(ROWSCA, M, RONE) 433 CALL SMUMPS_INITREAL(COLSCA, N, RONE) 434 ENDIF 435 ITDRPTR = 1 436 ITDCPTR = ITDRPTR + M 437C 438 ISRRPTR = ITDCPTR + N 439 OSRRPTR = ISRRPTR + IRSNDRCVVOL 440C 441 ISRCPTR = OSRRPTR + ORSNDRCVVOL 442 OSRCPTR = ISRCPTR + ICSNDRCVVOL 443C To avoid bound check errors... 444 IF(NUMPROCS == 1)THEN 445 OSRCPTR = OSRCPTR - 1 446 ISRCPTR = ISRCPTR - 1 447 OSRRPTR = OSRRPTR - 1 448 ISRRPTR = ISRRPTR - 1 449 ELSE 450 IF(IRSNDRCVVOL == 0) ISRRPTR = ISRRPTR - 1 451 IF(ORSNDRCVVOL == 0) OSRRPTR = OSRRPTR - 1 452 IF(ICSNDRCVVOL == 0) ISRCPTR = ISRCPTR - 1 453 IF(OCSNDRCVVOL == 0) OSRCPTR = OSRCPTR - 1 454 ENDIF 455C Check done outside 456C IF( OSRCPTR + OCSNDRCVVOL - 1 > RESZ) THEN 457C write(6,*) "Bora: NOTE: ", 458C & OSRCPTR + OCSNDRCVVOL - 1 , RESZ 459C write(6,*) "Bora: TODO. scimscaent_3 REPORT ERROR" 460C CALL flush(6) 461C ENDIF 462C computation starts 463 ITER = 1 464 DO WHILE (ITER.LE.NB1+NB2+NB3) 465C CLEAR temporary Dr and Dc 466 IF(NUMPROCS > 1) THEN 467 CALL SMUMPS_ZEROOUT(WRKRC(ITDRPTR),M, 468 & IWRK(IMYRPTR),INUMMYR) 469 CALL SMUMPS_ZEROOUT(WRKRC(ITDCPTR),N, 470 & IWRK(IMYCPTR),INUMMYC) 471 ELSE 472 CALL SMUMPS_INITREAL(WRKRC(ITDRPTR),M, RZERO) 473 CALL SMUMPS_INITREAL(WRKRC(ITDCPTR),N, RZERO) 474 ENDIF 475 IF((ITER.LE.NB1).OR.(ITER > NB1+NB2)) THEN 476C INF-NORM ITERATION 477 IF((ITER.EQ.1).OR.(OORANGEIND.EQ.1)) THEN 478 DO NZIND=1_8,NZ_loc 479 IR = IRN_loc(NZIND) 480 IC = JCN_loc(NZIND) 481 IF((IR.GE.1).AND.(IR.LE.M).AND. 482 & (IC.GE.1).AND.(IC.LE.N)) THEN 483 ELM = abs(A_loc(NZIND))*ROWSCA(IR)*COLSCA(IC) 484 IF(WRKRC(ITDRPTR-1+IR)<ELM) THEN 485 WRKRC(ITDRPTR-1+IR)= ELM 486 ENDIF 487 IF(WRKRC(ITDCPTR-1+IC)<ELM) THEN 488 WRKRC(ITDCPTR-1+IC)= ELM 489 ENDIF 490 ELSE 491 OORANGEIND = 1 492 ENDIF 493 ENDDO 494 ELSEIF(OORANGEIND.EQ.0) THEN 495 DO NZIND=1,NZ_loc 496 IR = IRN_loc(NZIND) 497 IC = JCN_loc(NZIND) 498 ELM = abs(A_loc(NZIND))*ROWSCA(IR)*COLSCA(IC) 499 IF(WRKRC(ITDRPTR-1+IR)<ELM) THEN 500 WRKRC(ITDRPTR-1+IR)= ELM 501 ENDIF 502 IF(WRKRC(ITDCPTR-1+IC)<ELM) THEN 503 WRKRC(ITDCPTR-1+IC)= ELM 504 ENDIF 505 ENDDO 506 ENDIF 507 IF(NUMPROCS > 1) THEN 508 CALL SMUMPS_DOCOMMINF(MYID, NUMPROCS, 509 & WRKRC(ITDCPTR), N, TAG_ITERS+ITER, 510 & ICSNDRCVNUM,IWRK(ICNGHBPRCS), 511 & ICSNDRCVVOL,IWRK(ICSNDRCVIA), IWRK(ICSNDRCVJA), 512 & WRKRC(ISRCPTR), 513 & OCSNDRCVNUM,IWRK(OCNGHBPRCS), 514 & OCSNDRCVVOL,IWRK(OCSNDRCVIA), IWRK(OCSNDRCVJA), 515 & WRKRC( OSRCPTR), 516 & IWRK(ISTATUS),IWRK(REQUESTS), 517 & COMM) 518C 519 CALL SMUMPS_DOCOMMINF(MYID, NUMPROCS, 520 & WRKRC(ITDRPTR), M, TAG_ITERS+2+ITER, 521 & IRSNDRCVNUM,IWRK(IRNGHBPRCS), 522 & IRSNDRCVVOL,IWRK(IRSNDRCVIA), IWRK(IRSNDRCVJA), 523 & WRKRC(ISRRPTR), 524 & ORSNDRCVNUM,IWRK(ORNGHBPRCS), 525 & ORSNDRCVVOL,IWRK(ORSNDRCVIA), IWRK(ORSNDRCVJA), 526 & WRKRC( OSRRPTR), 527 & IWRK(ISTATUS),IWRK(REQUESTS), 528 & COMM) 529 IF((EPS .GT. RZERO) .OR. 530 & (ITER.EQ.NB1).OR. 531 & ((ITER.EQ.NB1+NB2+NB3).AND. 532 & (NB1+NB3.GT.0))) THEN 533 INFERRROW = SMUMPS_ERRSCALOC(ROWSCA, 534 & WRKRC(ITDRPTR), M, 535 & IWRK(IMYRPTR),INUMMYR) 536C find error for the cols 537 INFERRCOL = SMUMPS_ERRSCALOC(COLSCA, 538 & WRKRC(ITDCPTR), N, 539 & IWRK(IMYCPTR),INUMMYC) 540C get max of those two errors 541 INFERRL = INFERRCOL 542 IF(INFERRROW > INFERRL ) THEN 543 INFERRL = INFERRROW 544 ENDIF 545C 546 CALL MPI_ALLREDUCE(INFERRL, INFERRG, 547 & 1, MPI_REAL, 548 & MPI_MAX, COMM, IERROR) 549 IF(INFERRG.LE.EPS) THEN 550 CALL SMUMPS_UPDATESCALE(COLSCA, 551 & WRKRC(ITDCPTR),N, 552 & IWRK(IMYCPTR),INUMMYC) 553 CALL SMUMPS_UPDATESCALE(ROWSCA, 554 & WRKRC(ITDRPTR),M, 555 & IWRK(IMYRPTR),INUMMYR) 556 IF(ITER .LE. NB1) THEN 557 ITER = NB1+1 558 CYCLE 559 ELSE 560 EXIT 561 ENDIF 562 ENDIF 563 ENDIF 564 ELSE 565C SINGLE PROCESSOR CASE: INF-NORM ERROR COMPUTATION 566 IF((EPS .GT. RZERO) .OR. 567 & (ITER.EQ.NB1).OR. 568 & ((ITER.EQ.NB1+NB2+NB3).AND. 569 & (NB1+NB3.GT.0))) THEN 570 INFERRROW = SMUMPS_ERRSCA1(ROWSCA, 571 & WRKRC(ITDRPTR), M) 572C find error for the cols 573 INFERRCOL = SMUMPS_ERRSCA1(COLSCA, 574 & WRKRC(ITDCPTR), N) 575C get max of those two errors 576 INFERRL = INFERRCOL 577 IF(INFERRROW > INFERRL) THEN 578 INFERRL = INFERRROW 579 ENDIF 580 INFERRG = INFERRL 581 IF(INFERRG.LE.EPS) THEN 582 CALL SMUMPS_UPSCALE1(COLSCA, WRKRC(ITDCPTR), N) 583 CALL SMUMPS_UPSCALE1(ROWSCA, WRKRC(ITDRPTR), M) 584 IF(ITER .LE. NB1) THEN 585 ITER = NB1+1 586 CYCLE 587 ELSE 588 EXIT 589 ENDIF 590 ENDIF 591 ENDIF 592 ENDIF 593 ELSE 594C WE HAVE ITER.GT.NB1 AND ITER.LE.NB1+NB2. 595C ONE-NORM ITERATION 596 IF((ITER .EQ.1).OR.(OORANGEIND.EQ.1))THEN 597 DO NZIND=1_8,NZ_loc 598 IR = IRN_loc(NZIND) 599 IC = JCN_loc(NZIND) 600 IF((IR.GE.1).AND.(IR.LE.M).AND. 601 & (IC.GE.1).AND.(IC.LE.N)) THEN 602 ELM = abs(A_loc(NZIND))*ROWSCA(IR)*COLSCA(IC) 603 WRKRC(ITDRPTR-1+IR) = WRKRC(ITDRPTR-1+IR) + ELM 604 WRKRC(ITDCPTR-1+IC) = WRKRC(ITDCPTR-1+IC) + ELM 605 ELSE 606 OORANGEIND = 1 607 ENDIF 608 ENDDO 609 ELSEIF(OORANGEIND.EQ.0) THEN 610 DO NZIND=1,NZ_loc 611 IR = IRN_loc(NZIND) 612 IC = JCN_loc(NZIND) 613 ELM = abs(A_loc(NZIND))*ROWSCA(IR)*COLSCA(IC) 614 WRKRC(ITDRPTR-1+IR) = WRKRC(ITDRPTR-1+IR) + ELM 615 WRKRC(ITDCPTR-1+IC) = WRKRC(ITDCPTR-1+IC) + ELM 616 ENDDO 617 ENDIF 618 IF(NUMPROCS > 1) THEN 619 CALL SMUMPS_DOCOMM1N(MYID, NUMPROCS, 620 & WRKRC(ITDCPTR), N, TAG_ITERS+ITER, 621 & ICSNDRCVNUM, IWRK(ICNGHBPRCS), 622 & ICSNDRCVVOL, IWRK(ICSNDRCVIA), IWRK(ICSNDRCVJA), 623 & WRKRC(ISRCPTR), 624 & OCSNDRCVNUM, IWRK(OCNGHBPRCS), 625 & OCSNDRCVVOL, IWRK(OCSNDRCVIA), IWRK(OCSNDRCVJA), 626 & WRKRC( OSRCPTR), 627 & IWRK(ISTATUS), IWRK(REQUESTS), 628 & COMM) 629C 630 CALL SMUMPS_DOCOMM1N(MYID, NUMPROCS, 631 & WRKRC(ITDRPTR), M, TAG_ITERS+2+ITER, 632 & IRSNDRCVNUM, IWRK(IRNGHBPRCS), 633 & IRSNDRCVVOL, IWRK(IRSNDRCVIA), IWRK(IRSNDRCVJA), 634 & WRKRC(ISRRPTR), 635 & ORSNDRCVNUM, IWRK(ORNGHBPRCS), 636 & ORSNDRCVVOL, IWRK(ORSNDRCVIA), IWRK(ORSNDRCVJA), 637 & WRKRC( OSRRPTR), 638 & IWRK(ISTATUS), IWRK(REQUESTS), 639 & COMM) 640 IF((EPS .GT. RZERO) .OR. 641 & ((ITER.EQ.NB1+NB2).AND. 642 & (NB2.GT.0))) THEN 643 ONEERRROW = SMUMPS_ERRSCALOC(ROWSCA, 644 & WRKRC(ITDRPTR), M, 645 & IWRK(IMYRPTR),INUMMYR) 646C find error for the cols 647 ONEERRCOL = SMUMPS_ERRSCALOC(COLSCA, 648 & WRKRC(ITDCPTR), N, 649 & IWRK(IMYCPTR),INUMMYC) 650C get max of those two errors 651 ONEERRL = ONEERRCOL 652 IF(ONEERRROW > ONEERRL ) THEN 653 ONEERRL = ONEERRROW 654 ENDIF 655C 656 CALL MPI_ALLREDUCE(ONEERRL, ONEERRG, 657 & 1, MPI_REAL, 658 & MPI_MAX, COMM, IERROR) 659 IF(ONEERRG.LE.EPS) THEN 660 CALL SMUMPS_UPDATESCALE(COLSCA, 661 & WRKRC(ITDCPTR),N, 662 & IWRK(IMYCPTR),INUMMYC) 663 CALL SMUMPS_UPDATESCALE(ROWSCA, 664 & WRKRC(ITDRPTR),M, 665 & IWRK(IMYRPTR),INUMMYR) 666 ITER = NB1+NB2+1 667 CYCLE 668 ENDIF 669 ENDIF 670 ELSE 671C SINGLE-PROCESSOR CASE: ONE-NORM ERROR COMPUTATION 672 IF((EPS .GT. RZERO) .OR. 673 & ((ITER.EQ.NB1+NB2).AND. 674 & (NB2.GT.0))) THEN 675 ONEERRROW = SMUMPS_ERRSCA1(ROWSCA, 676 & WRKRC(ITDRPTR), M) 677C find error for the cols 678 ONEERRCOL = SMUMPS_ERRSCA1(COLSCA, 679 & WRKRC(ITDCPTR), N) 680C get max of those two errors 681 ONEERRL = ONEERRCOL 682 IF(ONEERRROW > ONEERRL) THEN 683 ONEERRL = ONEERRROW 684 ENDIF 685 ONEERRG = ONEERRL 686 IF(ONEERRG.LE.EPS) THEN 687 CALL SMUMPS_UPSCALE1(COLSCA, WRKRC(ITDCPTR), N) 688 CALL SMUMPS_UPSCALE1(ROWSCA, WRKRC(ITDRPTR), M) 689 ITER = NB1+NB2+1 690 CYCLE 691 ENDIF 692 ENDIF 693 ENDIF 694 ENDIF 695 IF(NUMPROCS > 1) THEN 696 CALL SMUMPS_UPDATESCALE(COLSCA, WRKRC(ITDCPTR), N, 697 & IWRK(IMYCPTR),INUMMYC) 698 CALL SMUMPS_UPDATESCALE(ROWSCA, WRKRC(ITDRPTR), M, 699 & IWRK(IMYRPTR),INUMMYR) 700C 701 ELSE 702C SINGLE PROCESSOR CASE: Conv check and update of sca arrays 703 CALL SMUMPS_UPSCALE1(COLSCA, WRKRC(ITDCPTR), N) 704 CALL SMUMPS_UPSCALE1(ROWSCA, WRKRC(ITDRPTR), M) 705 ENDIF 706 ITER = ITER + 1 707 ENDDO 708 ONENORMERR = ONEERRG 709 INFNORMERR = INFERRG 710 IF(NUMPROCS > 1) THEN 711 CALL MPI_REDUCE(ROWSCA, WRKRC(1), M, MPI_REAL, 712 & MPI_MAX, 0, 713 & COMM, IERROR) 714 IF(MYID.EQ.0) THEN 715 DO I=1, M 716 ROWSCA(I) = WRKRC(I) 717 ENDDO 718 ENDIF 719C Scaling factors are printed 720C WRITE (6,*) MYID, 'ROWSCA=',ROWSCA 721C WRITE (6,*) MYID, 'COLSCA=',COLSCA 722C CALL FLUSH(6) 723c REduce the whole scaling factors to processor 0 of COMM 724 CALL MPI_REDUCE(COLSCA, WRKRC(1+M), N, MPI_REAL, 725 & MPI_MAX, 0, 726 & COMM, IERROR) 727 If(MYID.EQ.0) THEN 728 DO I=1, N 729 COLSCA(I) = WRKRC(I+M) 730 ENDDO 731 ENDIF 732 ENDIF 733 ENDIF 734 RETURN 735 END SUBROUTINE SMUMPS_SIMSCALEABSUNS 736C 737C 738C SEPARATOR: Another function begins 739C 740C 741 SUBROUTINE SMUMPS_SIMSCALEABSSYM(IRN_loc, JCN_loc, A_loc, NZ_loc, 742 & N, NUMPROCS, MYID, COMM, 743 & PARTVEC, 744 & RSNDRCVSZ, 745 & REGISTRE, 746 & IWRK, IWRKSZ, 747 & INTSZ, RESZ, OP, 748 & SCA, WRKRC, ISZWRKRC, 749 & NB1, NB2, NB3, EPS, 750 & ONENORMERR, INFNORMERR) 751C---------------------------------------------------------------------- 752C Input parameters: 753C N: size of matrix (sym matrix, square). 754C NUMPROCS, MYID, COMM: guess what are those 755C PARTVEC: row/col partvec to be filled when OP=1 756C RSNDRCVSZ:send recv sizes for row/col operations. 757C to be filled when OP=1 758C REGISTRE: to store some pointers (size etc). Its size is 12, 759C but we do not use all in this routine. 760C IWRK: working space. when OP=1 IWRKSZ.GE.2*MAXMN 761C when OP=2 INTSZ portion is used. Donc, IWRKSZ>INTSZ 762C when OP=2 763C IWRKSZ: size 764C INTSZ: to be computed when OP=1, necessary integer space to run 765C scaling algo when OP=2 766C RESZ: to be computed when OP=1, necessary real space to run 767C scaling algo when OP=2 768C OP: 769C =1 estimation of memory and construction of partvecs 770C writes into PARTVEC,RSNDRCVSZ,REGISTRE 771C does not access WRKRC, uses IWRK as workspace 772C computes INTSZ and RESZ. 773C =2 Compute scalings 774C restores pointers from REGISTRE, 775C stores communication structure in IWRK (from the start). 776C 777C SCA: space for row/col scaling factor; has size M 778C WRKRC: real working space. when OP=1, is not accessed. Donc, it 779C can be declared to be of size 1 at OP=1 call. 780C ISZWRKRC: size 781C SYM: is matrix symmetric 782C NB1, NB2, NB3: algo runs 783C NB1 iters of inf-norm (default 1/1), 784C NB2 iters of 1-norm (default 3/10), 785C NB3 iters of inf-norm (default 3/10). 786C in succession. 787C EPS: tolerance for concergence. 788C IF EPS < 0.R0 then does not test convergence. 789C See comments for the uns case above. 790C ONENORMERR : error in one norm scaling (see comments for the 791C uns case above), 792C INFNORMERR : error in inf norm scaling (see comments for the 793C uns case above). 794C--------------------------------------------------------------------- 795C On input: 796C OP=1==>Requirements 797C IWRKSZ.GE.2*MAXMN XXXX compare with uns variant. 798C PARTVEC of size N 799C SNDRCVSZ of size 2*NUMPROCS 800C REGISTRE of size 12 801C 802C OP=2==>Requirements 803C INTSZ .GE. REGISTRE(11) 804C RESZ .GE. REGISTRE(12) 805C--------------------------------------------------------------------- 806C On output: 807C SCA 808C at processor 0 of COMM: complete factors. 809C at other processors : only the SCA(i) and SCA(j) 810C for which there is a nonzero a_ij. 811C ONENORMERR : error in one norm scaling 812C = -1.0 if iter2=0. 813C INFNORMERR : error in inf norm scaling 814C = inf norm error at iter3 if iter3 > 0 815C = inf norm error at iter1 if iter1 > 0, iter3=0 816C = -1.0 if iter1=iter3=0 817C --------------------------------------------------------------------- 818C NOTE: some variables are named in such a way that they correspond 819C to the row variables in unsym case. They are used for both 820C row and col communications. 821C --------------------------------------------------------------------- 822C References: 823C The scaling algorithms are based on those discussed in 824C [1] D. Ruiz, "A scaling algorithm to equilibrate both rows and 825C columns norms in matrices", Tech. Rep. Rutherford 826C Appleton Laboratory, Oxon, UK and ENSEEIHT-IRIT, 827C Toulouse, France, RAL-TR-2001-034 and RT/APO/01/4, 2001. 828C [2] D. Ruiz and B. Ucar, "A symmetry preserving algorithm for 829C matrix scaling", in preparation as of Jan'08. 830C 831C The parallelization approach is based on discussion in 832C [3] P. R. Amestoy, I. S. Duff, D. Ruiz, and B. Ucar, "A parallel 833C matrix scaling algorithm", accepted for publication, 834C In proceedings of VECPAR'08-International Meeting-High 835C Performance Computing for Computational Science, Jan'08. 836C and was supported by ANR-SOLSTICE project (ANR-06-CIS6-010) 837C --------------------------------------------------------------------- 838 IMPLICIT NONE 839 INCLUDE 'mpif.h' 840 INTEGER(8) :: NZ_loc 841 INTEGER N, IWRKSZ, OP 842 INTEGER NUMPROCS, MYID, COMM 843 INTEGER INTSZ, RESZ 844 INTEGER IRN_loc(NZ_loc) 845 INTEGER JCN_loc(NZ_loc) 846 REAL A_loc(NZ_loc) 847 INTEGER PARTVEC(N), RSNDRCVSZ(2*NUMPROCS) 848 INTEGER IWRK(IWRKSZ) 849 INTEGER REGISTRE(12) 850 REAL SCA(N) 851 INTEGER ISZWRKRC 852 REAL WRKRC(ISZWRKRC) 853C LOCALS 854 INTEGER IRSNDRCVNUM, ORSNDRCVNUM 855 INTEGER IRSNDRCVVOL, ORSNDRCVVOL 856 INTEGER INUMMYR 857C IMPORTANT POINTERS 858 INTEGER IMYRPTR,IMYCPTR 859 INTEGER IRNGHBPRCS, IRSNDRCVIA,IRSNDRCVJA 860 INTEGER ORNGHBPRCS, ORSNDRCVIA,ORSNDRCVJA 861 INTEGER ISTATUS, REQUESTS, TMPWORK 862 INTEGER ITDRPTR, ISRRPTR, OSRRPTR 863 REAL ONENORMERR,INFNORMERR 864C FOR the scaling phase 865 INTEGER NB1, NB2, NB3 866 REAL EPS 867C Iteration vars 868 INTEGER ITER, IR, IC 869 INTEGER(8) :: NZIND 870 REAL ELM 871C COMM TAGS.... 872 INTEGER TAG_COMM_ROW 873 PARAMETER(TAG_COMM_ROW=101) 874 INTEGER TAG_ITERS 875 PARAMETER(TAG_ITERS=102) 876C FUNCTIONS 877 EXTERNAL SMUMPS_CREATEPARTVECSYM, 878 & SMUMPS_NUMVOLSNDRCVSYM, 879 & SMUMPS_SETUPCOMMSSYM, 880 & SMUMPS_FINDNUMMYROWCOLSYM, 881 & SMUMPS_CHKCONVGLOSYM, 882 & SMUMPS_CHK1CONV, 883 & SMUMPS_FILLMYROWCOLINDICESSYM, 884 & SMUMPS_DOCOMMINF, 885 & SMUMPS_DOCOMM1N, 886 & SMUMPS_INITREAL, 887 & SMUMPS_INITREALLST 888 INTEGER SMUMPS_CHKCONVGLOSYM 889 INTEGER SMUMPS_CHK1CONV 890 REAL SMUMPS_ERRSCALOC 891 REAL SMUMPS_ERRSCA1 892 INTRINSIC abs 893 REAL RONE, RZERO 894 PARAMETER(RONE=1.0E0,RZERO=0.0E0) 895C TMP VARS 896 INTEGER INTSZR 897 INTEGER MAXMN 898 INTEGER I, IERROR 899 REAL ONEERRL, ONEERRG 900 REAL INFERRL, INFERRG 901 INTEGER OORANGEIND 902 OORANGEIND = 0 903 INFERRG = -RONE 904 ONEERRG = -RONE 905 MAXMN = N 906 IF(OP == 1) THEN 907 IF(NUMPROCS > 1) THEN 908C Check done outside 909C IF(IWRKSZ.LT.2*MAXMN) THEN ERROR.... 910 CALL SMUMPS_CREATEPARTVECSYM(MYID, NUMPROCS, COMM, 911 & IRN_loc, JCN_loc, NZ_loc, 912 & PARTVEC, N, 913 & IWRK, IWRKSZ) 914C 915 CALL SMUMPS_NUMVOLSNDRCVSYM(MYID, NUMPROCS, N, PARTVEC, 916 & NZ_loc, IRN_loc, JCN_loc, IRSNDRCVNUM,IRSNDRCVVOL, 917 & ORSNDRCVNUM, ORSNDRCVVOL, 918 & IWRK,IWRKSZ, 919 & RSNDRCVSZ(1), RSNDRCVSZ(1+NUMPROCS), COMM) 920C 921 CALL SMUMPS_FINDNUMMYROWCOLSYM(MYID, NUMPROCS, COMM, 922 & IRN_loc, JCN_loc, NZ_loc, 923 & PARTVEC, N, 924 & INUMMYR, 925 & IWRK, IWRKSZ) 926C 927 INTSZR = IRSNDRCVNUM + ORSNDRCVNUM + 928 & IRSNDRCVVOL + ORSNDRCVVOL + 929 & 2*(NUMPROCS+1) + INUMMYR 930 INTSZ = INTSZR + N + 931 & (MPI_STATUS_SIZE +1) * NUMPROCS 932 ELSE 933C NUMPROCS IS 1 934 IRSNDRCVNUM = 0 935 ORSNDRCVNUM = 0 936 IRSNDRCVVOL = 0 937 ORSNDRCVVOL = 0 938 INUMMYR = 0 939 INTSZ = 0 940 ENDIF 941C CALCULATE NECESSARY REAL SPACE 942 RESZ = N + IRSNDRCVVOL + ORSNDRCVVOL 943C write(6,*) 'Bora :', RESZ, N, IRSNDRCVVOL, ORSNDRCVVOL 944C CALL flush(6) 945C CALCULATE NECESSARY INT SPACE 946C The last maxmn is tmpwork for setup comm and fillmyrowcol 947 REGISTRE(1) = IRSNDRCVNUM 948 REGISTRE(2) = ORSNDRCVNUM 949 REGISTRE(3) = IRSNDRCVVOL 950 REGISTRE(4) = ORSNDRCVVOL 951 REGISTRE(9) = INUMMYR 952 REGISTRE(11) = INTSZ 953 REGISTRE(12) = RESZ 954 ELSE 955C else of op=1. That is op=2 now. 956C restore the numbers 957 IRSNDRCVNUM = REGISTRE(1) 958 ORSNDRCVNUM = REGISTRE(2) 959 IRSNDRCVVOL = REGISTRE(3) 960 ORSNDRCVVOL = REGISTRE(4) 961 INUMMYR = REGISTRE(9) 962 IF(NUMPROCS > 1) THEN 963C Check done outsize 964C IF(INTSZ < REGISTRE(11)) THEN ERROR 965C IF(RESZ < REGISTRE(12)) THEN ERROR 966C Fill up myrows and my colsX 967 CALL SMUMPS_FILLMYROWCOLINDICESSYM(MYID, NUMPROCS,COMM, 968 & IRN_loc, JCN_loc, NZ_loc, 969 & PARTVEC, N, 970 & IWRK(1), INUMMYR, 971 & IWRK(1+INUMMYR), IWRKSZ-INUMMYR) 972 IMYRPTR = 1 973 IMYCPTR = IMYRPTR + INUMMYR 974C Set up comm and run. 975C set pointers in iwrk (3 parts) 976C 977C ROWS [---------------------------------------------] 978 IRNGHBPRCS = IMYCPTR 979 IRSNDRCVIA = IRNGHBPRCS+IRSNDRCVNUM 980 IRSNDRCVJA = IRSNDRCVIA + NUMPROCS+1 981 ORNGHBPRCS = IRSNDRCVJA + IRSNDRCVVOL 982 ORSNDRCVIA = ORNGHBPRCS + ORSNDRCVNUM 983 ORSNDRCVJA = ORSNDRCVIA + NUMPROCS + 1 984C MPI [-----------------] 985 REQUESTS = ORSNDRCVJA + ORSNDRCVVOL 986 ISTATUS = REQUESTS + NUMPROCS 987C TMPWRK [-----------------] 988 TMPWORK = ISTATUS + MPI_STATUS_SIZE * NUMPROCS 989 CALL SMUMPS_SETUPCOMMSSYM(MYID, NUMPROCS, N, PARTVEC, 990 & NZ_loc, IRN_loc, JCN_loc, 991 & IRSNDRCVNUM, IRSNDRCVVOL, 992 & IWRK(IRNGHBPRCS),IWRK(IRSNDRCVIA),IWRK(IRSNDRCVJA), 993 & ORSNDRCVNUM, ORSNDRCVVOL, 994 & IWRK(ORNGHBPRCS),IWRK(ORSNDRCVIA),IWRK(ORSNDRCVJA), 995 & RSNDRCVSZ(1), RSNDRCVSZ(1+NUMPROCS), 996 & IWRK(TMPWORK), 997 & IWRK(ISTATUS), IWRK(REQUESTS), 998 & TAG_COMM_ROW, COMM) 999 CALL SMUMPS_INITREAL(SCA, N, RZERO) 1000 CALL SMUMPS_INITREALLST(SCA, N, 1001 & IWRK(IMYRPTR),INUMMYR, RONE) 1002 ELSE 1003 CALL SMUMPS_INITREAL(SCA, N, RONE) 1004 ENDIF 1005 ITDRPTR = 1 1006 ISRRPTR = ITDRPTR + N 1007 OSRRPTR = ISRRPTR + IRSNDRCVVOL 1008C 1009C To avoid bound check errors... 1010 IF(NUMPROCS == 1)THEN 1011 OSRRPTR = OSRRPTR - 1 1012 ISRRPTR = ISRRPTR - 1 1013 ELSE 1014 IF(IRSNDRCVVOL == 0) ISRRPTR = ISRRPTR - 1 1015 IF(ORSNDRCVVOL == 0) OSRRPTR = OSRRPTR - 1 1016 ENDIF 1017C computation starts 1018 ITER = 1 1019 DO WHILE(ITER.LE.NB1+NB2+NB3) 1020C CLEAR temporary Dr and Dc 1021 IF(NUMPROCS > 1) THEN 1022 CALL SMUMPS_ZEROOUT(WRKRC(ITDRPTR),N, 1023 & IWRK(IMYRPTR),INUMMYR) 1024 ELSE 1025 CALL SMUMPS_INITREAL(WRKRC(ITDRPTR),N, RZERO) 1026 ENDIF 1027 IF((ITER.LE.NB1).OR.(ITER > NB1+NB2)) THEN 1028C INF-NORM ITERATION 1029 IF((ITER.EQ.1).OR.(OORANGEIND.EQ.1)) THEN 1030 DO NZIND=1_8,NZ_loc 1031 IR = IRN_loc(NZIND) 1032 IC = JCN_loc(NZIND) 1033 IF((IR.GE.1).AND.(IR.LE.N).AND. 1034 & (IC.GE.1).AND.(IC.LE.N)) THEN 1035 ELM = abs(A_loc(NZIND))*SCA(IR)*SCA(IC) 1036 IF(WRKRC(ITDRPTR-1+IR)<ELM) THEN 1037 WRKRC(ITDRPTR-1+IR)= ELM 1038 ENDIF 1039 IF(WRKRC(ITDRPTR-1+IC)<ELM) THEN 1040 WRKRC(ITDRPTR-1+IC)= ELM 1041 ENDIF 1042 ELSE 1043 OORANGEIND = 1 1044 ENDIF 1045 ENDDO 1046 ELSEIF(OORANGEIND.EQ.0) THEN 1047 DO NZIND=1_8,NZ_loc 1048 IR = IRN_loc(NZIND) 1049 IC = JCN_loc(NZIND) 1050 ELM = abs(A_loc(NZIND))*SCA(IR)*SCA(IC) 1051 IF(WRKRC(ITDRPTR-1+IR)<ELM) THEN 1052 WRKRC(ITDRPTR-1+IR)= ELM 1053 ENDIF 1054 IF(WRKRC(ITDRPTR-1+IC)<ELM) THEN 1055 WRKRC(ITDRPTR-1+IC)= ELM 1056 ENDIF 1057 ENDDO 1058 ENDIF 1059 IF(NUMPROCS > 1) THEN 1060 CALL SMUMPS_DOCOMMINF(MYID, NUMPROCS, 1061 & WRKRC(ITDRPTR), N, TAG_ITERS+2+ITER, 1062 & IRSNDRCVNUM,IWRK(IRNGHBPRCS), 1063 & IRSNDRCVVOL,IWRK(IRSNDRCVIA), IWRK(IRSNDRCVJA), 1064 & WRKRC(ISRRPTR), 1065 & ORSNDRCVNUM,IWRK(ORNGHBPRCS), 1066 & ORSNDRCVVOL,IWRK(ORSNDRCVIA), IWRK(ORSNDRCVJA), 1067 & WRKRC( OSRRPTR), 1068 & IWRK(ISTATUS),IWRK(REQUESTS), 1069 & COMM) 1070 IF((EPS .GT. RZERO) .OR. 1071 & (ITER.EQ.NB1).OR. 1072 & ((ITER.EQ.NB1+NB2+NB3).AND. 1073 & (NB1+NB3.GT.0))) THEN 1074 INFERRL = SMUMPS_ERRSCALOC(SCA, 1075 & WRKRC(ITDRPTR), N, 1076 & IWRK(IMYRPTR),INUMMYR) 1077 CALL MPI_ALLREDUCE(INFERRL, INFERRG, 1078 & 1, MPI_REAL, 1079 & MPI_MAX, COMM, IERROR) 1080 IF(INFERRG.LE.EPS) THEN 1081 CALL SMUMPS_UPDATESCALE(SCA, WRKRC(ITDRPTR), N, 1082 & IWRK(IMYRPTR),INUMMYR) 1083 IF(ITER .LE. NB1) THEN 1084 ITER = NB1+1 1085 CYCLE 1086 ELSE 1087 EXIT 1088 ENDIF 1089 ENDIF 1090 ENDIF 1091 ELSE 1092C SINGLE PROCESSOR CASE: INF-NORM ERROR COMPUTATION 1093 IF((EPS .GT. RZERO) .OR. 1094 & (ITER.EQ.NB1).OR. 1095 & ((ITER.EQ.NB1+NB2+NB3).AND. 1096 & (NB1+NB3.GT.0))) THEN 1097 INFERRL = SMUMPS_ERRSCA1(SCA, 1098 & WRKRC(ITDRPTR), N) 1099 INFERRG = INFERRL 1100 IF(INFERRG.LE.EPS) THEN 1101 CALL SMUMPS_UPSCALE1(SCA, WRKRC(ITDRPTR), N) 1102 IF(ITER .LE. NB1) THEN 1103 ITER = NB1+1 1104 CYCLE 1105 ELSE 1106 EXIT 1107 ENDIF 1108 ENDIF 1109 ENDIF 1110 ENDIF 1111 ELSE 1112C WE HAVE ITER.GT.NB1 AND ITER.LE.NB1+NB2. 1113C ONE-NORM ITERATION 1114 IF((ITER.EQ.1).OR.(OORANGEIND.EQ.1))THEN 1115 DO NZIND=1_8,NZ_loc 1116 IR = IRN_loc(NZIND) 1117 IC = JCN_loc(NZIND) 1118 IF((IR.GE.1).AND.(IR.LE.N).AND. 1119 & (IC.GE.1).AND.(IC.LE.N)) THEN 1120 ELM = abs(A_loc(NZIND))*SCA(IR)*SCA(IC) 1121 WRKRC(ITDRPTR-1+IR) = WRKRC(ITDRPTR-1+IR) + ELM 1122 IF(IR.NE.IC) THEN 1123 WRKRC(ITDRPTR-1+IC) = 1124 & WRKRC(ITDRPTR-1+IC) + ELM 1125 ENDIF 1126 ELSE 1127 OORANGEIND = 1 1128 ENDIF 1129 ENDDO 1130 ELSEIF(OORANGEIND.EQ.0)THEN 1131 DO NZIND=1_8,NZ_loc 1132 IR = IRN_loc(NZIND) 1133 IC = JCN_loc(NZIND) 1134 ELM = abs(A_loc(NZIND))*SCA(IR)*SCA(IC) 1135 WRKRC(ITDRPTR-1+IR) = WRKRC(ITDRPTR-1+IR) + ELM 1136 IF(IR.NE.IC) THEN 1137 WRKRC(ITDRPTR-1+IC) = WRKRC(ITDRPTR-1+IC) + ELM 1138 ENDIF 1139 ENDDO 1140 ENDIF 1141 IF(NUMPROCS > 1) THEN 1142 CALL SMUMPS_DOCOMM1N(MYID, NUMPROCS, 1143 & WRKRC(ITDRPTR), N, TAG_ITERS+2+ITER, 1144 & IRSNDRCVNUM, IWRK(IRNGHBPRCS), 1145 & IRSNDRCVVOL, IWRK(IRSNDRCVIA), IWRK(IRSNDRCVJA), 1146 & WRKRC(ISRRPTR), 1147 & ORSNDRCVNUM, IWRK(ORNGHBPRCS), 1148 & ORSNDRCVVOL, IWRK(ORSNDRCVIA), IWRK(ORSNDRCVJA), 1149 & WRKRC( OSRRPTR), 1150 & IWRK(ISTATUS), IWRK(REQUESTS), 1151 & COMM) 1152 IF((EPS .GT. RZERO) .OR. 1153 & ((ITER.EQ.NB1+NB2).AND. 1154 & (NB2.GT.0))) THEN 1155 ONEERRL = SMUMPS_ERRSCALOC(SCA, 1156 & WRKRC(ITDRPTR), N, 1157 & IWRK(IMYRPTR),INUMMYR) 1158C mpi allreduce. 1159 CALL MPI_ALLREDUCE(ONEERRL, ONEERRG, 1160 & 1, MPI_REAL, 1161 & MPI_MAX, COMM, IERROR) 1162 IF(ONEERRG.LE.EPS) THEN 1163 CALL SMUMPS_UPDATESCALE(SCA, WRKRC(ITDRPTR), N, 1164 & IWRK(IMYRPTR),INUMMYR) 1165 ITER = NB1+NB2+1 1166 CYCLE 1167 ENDIF 1168 ENDIF 1169 ELSE 1170C SINGLE-PROCESSOR CASE: ONE-NORM ERROR COMPUTATION 1171 IF((EPS .GT. RZERO) .OR. 1172 & ((ITER.EQ.NB1+NB2).AND. 1173 & (NB2.GT.0))) THEN 1174 ONEERRL = SMUMPS_ERRSCA1(SCA, 1175 & WRKRC(ITDRPTR), N) 1176 ONEERRG = ONEERRL 1177 IF(ONEERRG.LE.EPS) THEN 1178 CALL SMUMPS_UPSCALE1(SCA, WRKRC(ITDRPTR), N) 1179 ITER = NB1+NB2+1 1180 CYCLE 1181 ENDIF 1182 ENDIF 1183 ENDIF 1184 ENDIF 1185 IF(NUMPROCS > 1) THEN 1186 CALL SMUMPS_UPDATESCALE(SCA, WRKRC(ITDRPTR), N, 1187 & IWRK(IMYRPTR),INUMMYR) 1188 ELSE 1189 CALL SMUMPS_UPSCALE1(SCA, WRKRC(ITDRPTR), N) 1190 ENDIF 1191 ITER = ITER + 1 1192 ENDDO 1193 ONENORMERR = ONEERRG 1194 INFNORMERR = INFERRG 1195 IF(NUMPROCS > 1) THEN 1196 CALL MPI_REDUCE(SCA, WRKRC(1), N, MPI_REAL, 1197 & MPI_MAX, 0, 1198 & COMM, IERROR) 1199 IF(MYID.EQ.0) THEN 1200 DO I=1, N 1201 SCA(I) = WRKRC(I) 1202 ENDDO 1203 ENDIF 1204 ENDIF 1205 ENDIF 1206 RETURN 1207 END SUBROUTINE SMUMPS_SIMSCALEABSSYM 1208