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