V3FIT
ptrd_mod.f90
1  module ptrd_mod
2  use stel_kinds
3  USE v3_utilities, ONLY: assert
4  implicit none
5 
6 #if defined(MPI_OPT)
7  contains
8 
9  subroutine pdtrdf(m,nblock,Amat,Bmat,Cmat,ipivmat, desc)
10  implicit none
11 !
12 ! perform LU factorization of block tridiagonal system
13 !
14 ! [A1 C1 ] [ L1 ] [ I U1 ]
15 ! [B1 A2 C2 ] [ B1 L2 ] [ I U2 ]
16 ! [ B2 A3 C3 ] = [ B2 L3 ] [ I U3 ]
17 ! [ B4 A4 ] [ B3 L4 ] [ I ]
18 !
19 !
20 ! L1 = A1
21 ! L1*U1 = C1 => U1 = L1\C1
22 ! B1*U1 + L2 = A2 => L2 = A2 - B1*U1
23 ! L2*U2 = C2 => U2 = L2\C2
24 !
25 ! Lk = Ak - B_{k-1}*U_{k-1}
26 ! Uk = Lk\Ck
27 !
28 !
29  REAL(rprec), parameter :: one = 1.0d0
30  REAL(rprec), parameter :: zero = 0.0d0
31 
32  integer, intent(in) :: m, nblock
33  REAL(rprec), target, dimension(:,:), intent(inout) :: Amat
34  REAL(rprec), target, dimension(:,:), intent(inout) :: Bmat
35  REAL(rprec), target, dimension(:,:), intent(inout) :: Cmat
36  integer, dimension(:,:), target, intent(inout) :: ipivmat
37  integer, dimension(:), intent(in) :: desc
38 
39 
40  integer :: k, ia,ja, ib,jb, iu,ju, ic,jc
41  integer :: info, nrhs, mm,nn
42  integer, dimension(:), pointer :: ipiv
43  REAL(rprec) :: alpha, beta
44  REAL(rprec), dimension(:), pointer :: Ak, Bkm1, Ukm1
45  REAL(rprec), dimension(:), pointer :: Lk, Ck
46 
47 
48  nullify( ak )
49  nullify( bkm1 )
50  nullify( ukm1 )
51  nullify( lk )
52  nullify( ck )
53 
54 
55  do k=1,nblock
56  if (k .ge. 2) then
57 ! -------------------------
58 ! Ak = Ak - B_{k-1}*U_{k-1}
59 ! -------------------------
60  ak => amat(:,k)
61  bkm1 => bmat(:,k-1)
62  ukm1 => cmat(:,k-1)
63  alpha = -one
64  beta = one
65  ib = 1
66  jb = 1
67  iu = 1
68  ju = 1
69  ia = 1
70  ja = 1
71 
72  call pdgemm( 'N', 'N', m,m,m, &
73  & alpha, bkm1,ib,jb,desc, &
74  & ukm1,iu,ju,desc, &
75  & beta, ak, ia,ja, desc )
76  endif
77 
78 ! ------------
79 ! Uk = Lk \ Ck
80 ! ------------
81  lk => amat(:,k)
82  ipiv => ipivmat(:,k)
83  ia = 1
84  ja = 1
85  info = 0
86  mm = m
87  nn = m
88  call pdgetrf(mm,nn,lk,ia,ja,desc,ipiv,info)
89  call assert(info.eq.0,'pdtrdf: pdgetrf return info != 0')
90 
91  if (k .le. (nblock-1)) then
92  nrhs = m
93  ic = 1
94  jc = 1
95  ck => cmat(:,k)
96  call pdgetrs('No',m,nrhs,lk,ia,ja,desc,ipiv, &
97  & ck,ic,jc,desc, info)
98  call assert(info.eq.0,'pdtrdf: pdgetrs return info != 0')
99  endif
100 
101 
102  enddo
103  return
104  end subroutine pdtrdf
105 
106  subroutine pdtrds(m,nblock,Amat,Bmat,Cmat,ipivmat,desc, &
107  & nrhs, Rrhs_in,ir,jr,descRrhs_in)
108  use descriptor_mod, ir1=>ir, jr1=>jr, mm1=>mm, info1=>info
109  implicit none
110 
111 
112 !
113 ! use LU factorization of block tridiagonal system
114 !
115 ! [A1 C1 ] [ L1 ] [ I U1 ]
116 ! [B1 A2 C2 ] [ B1 L2 ] [ I U2 ]
117 ! [ B2 A3 C3 ] = [ B2 L3 ] [ I U3 ]
118 ! [ B4 A4 ] [ B3 L4 ] [ I ]
119 !
120 !
121 !
122 !
123  integer, parameter :: idebug = 0
124  REAL(rprec), parameter :: one = 1.0d0
125  REAL(rprec), parameter :: zero = 0.0d0
126 
127  integer, intent(in) :: m, nblock
128  REAL(rprec), dimension(:,:), target, intent(inout) :: Amat
129  REAL(rprec), dimension(:,:), target, intent(inout) :: Bmat
130  REAL(rprec), dimension(:,:), target, intent(inout) :: Cmat
131  integer, dimension(:,:), target, intent(inout) :: ipivmat
132  integer, dimension(:), intent(in) :: desc
133 
134  integer, intent(in) :: nrhs
135  integer, intent(in) :: ir,jr
136  REAL(rprec), dimension(:) :: Rrhs_in
137  integer, dimension(:), intent(in) :: descRrhs_in
138 
139 
140  integer :: k, info, mm,nn,kk
141  integer :: ia,ja, irk, jrk, iy,jy, ix,jx, iu,ju, ib,jb
142  REAL(rprec) :: alpha, beta
143 
144  REAL(rprec), pointer, dimension(:) :: Lk,Bkm1,Uk
145  REAL(rprec), pointer, dimension(:) :: Rrhsk,Rrhskm1,Rrhskp1
146  integer, pointer, dimension(:) :: ipiv
147 
148 
149  integer :: inc1,inc2,iblock,irhs
150 
151  logical, parameter :: use_pdcopy = .false.
152  REAL(rprec), dimension(:,:), target, allocatable :: Rrhs
153  integer, dimension(DLEN_) :: descRrhs
154  INTEGER :: NUMROC
155  EXTERNAL :: numroc
156 
157 ! -----------------------------------------------------
158 ! change storage to nblock copies of m by nrhs matrices
159 ! also to avoid alignment constraints in scalapack
160 ! -----------------------------------------------------
161  icontxt = descrrhs_in(ctxt_)
162  mb = descrrhs_in(mb_)
163  nb = 1
164  rsrc = 0
165  csrc = 0
166 
167  call blacs_gridinfo( icontxt, nprow,npcol,myrow,mycol)
168  locp = numroc( m, mb, myrow,rsrc,nprow)
169  locq = numroc( nrhs,nb,mycol,csrc,npcol)
170  lld = max(1,locp)
171  call descinit(descrrhs,m,nrhs,mb,nb,rsrc,csrc,icontxt,lld,info)
172  call assert(info.eq.0,'pdtrds: descinit return info != 0')
173 
174  ineed = max(1,locp*locq)
175  allocate( rrhs(ineed,nblock),stat=ierr)
176  call assert(ierr.eq.0,'pdtrds: alloc Rrhs,ierr != 0')
177  rrhs = 0.0d0
178 
179  nullify( rrhsk )
180  nullify( rrhskm1 )
181  nullify( rrhskp1 )
182 
183 ! ---------
184 ! copy data
185 ! ---------
186 
187  if (use_pdcopy) then
188 
189  do iblock=1,nblock
190  do irhs=1,nrhs
191  ia = (ir-1) + 1 + (iblock-1)*m
192  ja = (jr-1) + irhs
193  inc1 = 1
194 
195  ib = 1
196  jb = irhs
197  inc2 = 1
198 
199  rrhsk => rrhs(:,iblock)
200  call pdcopy(m, rrhs_in,ia,ja,descrrhs_in,inc1, &
201  & rrhsk,ib,jb,descrrhs,inc2 )
202  enddo
203  enddo
204 
205  else
206 
207  do iblock=1,nblock
208  irhs = 1
209  ia = (ir-1) + 1 + (iblock-1)*m
210  ja = (jr-1) + irhs
211 
212  ib = 1
213  jb = irhs
214  rrhsk => rrhs(:,iblock)
215  alpha = 1.0d0
216  beta = 0.0d0
217  call pdgeadd( 'N',m,nrhs,alpha,rrhs_in,ia,ja,descrrhs_in, &
218  & beta, rrhsk,ib,jb,descrrhs )
219  enddo
220 
221  endif
222 
223 ! -----------------
224 ! L*U*x = r
225 ! (1) solve L*y = r
226 ! (2) solve U*x = y
227 ! -----------------
228  nullify(lk)
229  nullify(bkm1)
230  nullify(uk)
231  nullify(ipiv)
232 
233  icontxt = desc(ctxt_)
234  call blacs_gridinfo(icontxt, nprow,npcol,myrow,mycol)
235  isroot = (myrow.eq.0).and.(mycol.eq.0)
236 
237 ! --------------------------------
238 ! (1) solve L*y = r
239 !
240 ! [L1 ] [ y1 ] [ r1 ]
241 ! [B1 L2 ] [ y2 ] [ r2 ]
242 ! [ B2 L3 ] [ y3 ] = [ r3 ]
243 ! [ B3 L4 ] [ y4 ] [ r4 ]
244 !
245 !
246 ! y1 = L1\r1
247 ! y2 = L2\( r2 - B1*y1 )
248 ! y3 = L3\( r3 - B2*y2 )
249 ! y4 = L4\( r4 - B3*y3 )
250 !
251 ! yk = Lk\( rk - B_{k-1}*y_{k-1} )
252 ! --------------------------------
253  if (isroot .and. (idebug.ge.1)) then
254  write(*,*) 'pdtrds 77: m,nblock,nrhs ',m,nblock,nrhs
255  write(*,*) 'descRrhs_in(M_) ',descrrhs_in(m_)
256  write(*,*) 'descRrhs_in(N_) ',descrrhs_in(n_)
257  write(*,*) 'descRrhs_in(MB_) ',descrrhs_in(mb_)
258  write(*,*) 'descRrhs_in(NB_) ',descrrhs_in(nb_)
259  endif
260 
261  do k=1,nblock
262  if (k .ge. 2) then
263 ! --------------------------
264 ! rk <- rk - B_{k-1}*y_{k-1}
265 ! --------------------------
266  bkm1 => bmat(:,k-1)
267 
268  alpha = -one
269  beta = one
270  mm = m
271  nn = nrhs
272  kk = m
273  ib = 1
274  jb = 1
275 
276  iy = 1
277  jy = 1
278  irk = 1
279  jrk = 1
280  rrhsk => rrhs(:,k)
281  rrhskm1 => rrhs(:,k-1)
282 
283  call pdgemm( 'N', 'N', mm,nn,kk, &
284  & alpha, bkm1, ib,jb, desc, &
285  & rrhskm1, iy,jy, descrrhs, &
286  & beta, rrhsk, irk,jrk,descrrhs )
287  endif
288 
289 ! ------------
290 ! yk = Lk \ rk
291 ! ------------
292 
293  lk => amat(:,k)
294  ia = 1
295  ja = 1
296  ipiv => ipivmat(:,k)
297 
298  rrhsk => rrhs(:,k)
299  jrk = 1
300  irk = 1
301 
302  info = 0
303  if (isroot .and. (idebug.ge.1)) then
304  write(*,*) 'pdtrds 106: k,irk,jrk ',k,irk,jrk
305  endif
306  call pdgetrs( 'N', m,nrhs, lk,ia,ja,desc, ipiv, &
307  & rrhsk,irk,jrk,descrrhs,info)
308  call assert(info.eq.0,'pdtrds: pdgetrs return info != 0')
309 
310  enddo
311 
312 
313 ! (2) solve U*x = y
314 !
315 ! [I U1 ] [ x1 ] [ y1 ]
316 ! [ I U2 ] [ x2 ] [ y2 ]
317 ! [ I U3 ] [ x3 ] = [ y3 ]
318 ! [ I ] [ x4 ] [ y4 ]
319 !
320 !
321 ! x4 = y4
322 ! x3 = y3 - U3*y4
323 ! x2 = y2 - U2*y3
324 ! x1 = y1 - U1*y2
325 !
326 ! xk = yk - Uk*x_{k+1}
327 !
328 
329  do k=nblock-1,1,-1
330 ! --------------------
331 ! xk = yk - Uk*x_{k+1}
332 ! --------------------
333  alpha = -one
334  beta = one
335 
336  mm = m
337  nn = nrhs
338  kk = m
339  uk => cmat(:,k)
340  iu = 1
341  ju = 1
342 
343  ix = 1
344  jx = 1
345  irk = 1
346  jrk = 1
347 
348  rrhskp1 => rrhs(:,k+1)
349  rrhsk => rrhs(:,k)
350 
351  call pdgemm( 'N','N', mm,nn,kk, &
352  & alpha, uk,iu,ju,desc, &
353  & rrhskp1,ix,jx,descrrhs, &
354  & beta, rrhsk,irk,jrk,descrrhs )
355  enddo
356 
357 ! ----------------
358 ! copy results out
359 ! ----------------
360 
361  if (use_pdcopy) then
362 
363  do iblock=1,nblock
364  do irhs=1,nrhs
365  ia = (ir-1) + 1 + (iblock-1)*m
366  ja = (jr-1) + irhs
367  inc1 = 1
368 
369  ib = 1
370  jb = irhs
371  inc2 = 1
372 
373  rrhsk => rrhs(:,iblock)
374  call pdcopy(m,rrhsk,ib,jb,descrrhs,inc2, &
375  & rrhs_in,ia,ja,descrrhs_in,inc1)
376  enddo
377  enddo
378 
379  else
380 
381  do iblock=1,nblock
382  irhs = 1
383  ia = (ir-1) + 1 + (iblock-1)*m
384  ja = (jr-1) + irhs
385 
386  ib = 1
387  jb = irhs
388  rrhsk => rrhs(:,iblock)
389  alpha = 1.0d0
390  beta = 0.0d0
391 
392  call pdgeadd('N',m,nrhs,alpha,rrhsk,ib,jb,descrrhs, &
393  & beta,rrhs_in,ia,ja,descrrhs_in)
394  enddo
395 
396  endif
397  deallocate( rrhs, stat=ierr)
398  call assert(ierr.eq.0,'pdtrds: dealloc Rhs ')
399 
400 
401  return
402  end subroutine pdtrds
403 
404 #if defined(NEED_TOOLS)
405  SUBROUTINE pdgetrs( TRANS, N, NRHS, A, IA, JA, DESCA, IPIV, B, &
406  & IB, JB, DESCB, INFO )
407 !
408 ! -- ScaLAPACK routine (version 1.7) --
409 ! University of Tennessee, Knoxville, Oak Ridge National Laboratory,
410 ! and University of California, Berkeley.
411 ! May 1, 1997
412 !
413 ! .. Scalar Arguments ..
414  CHARACTER TRANS
415  INTEGER IA, IB, INFO, JA, JB, N, NRHS
416 ! ..
417 ! .. Array Arguments ..
418  INTEGER DESCA( * ), DESCB( * ), IPIV( * )
419  REAL(dp) A( * ), B( * )
420 ! ..
421 !
422 ! Purpose
423 ! =======
424 !
425 ! PDGETRS solves a system of distributed linear equations
426 !
427 ! op( sub( A ) ) * X = sub( B )
428 !
429 ! with a general N-by-N distributed matrix sub( A ) using the LU
430 ! factorization computed by PDGETRF.
431 ! sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1), op( A ) = A or A**T and
432 ! sub( B ) denotes B(IB:IB+N-1,JB:JB+NRHS-1).
433 !
434 ! Notes
435 ! =====
436 !
437 ! Each global data object is described by an associated description
438 ! vector. This vector stores the information required to establish
439 ! the mapping between an object element and its corresponding process
440 ! and memory location.
441 !
442 ! Let A be a generic term for any 2D block cyclicly distributed array.
443 ! Such a global array has an associated description vector DESCA.
444 ! In the following comments, the character _ should be read as
445 ! "of the global array".
446 !
447 ! NOTATION STORED IN EXPLANATION
448 ! --------------- -------------- --------------------------------------
449 ! DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
450 ! DTYPE_A = 1.
451 ! CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
452 ! the BLACS process grid A is distribu-
453 ! ted over. The context itself is glo-
454 ! bal, but the handle (the integer
455 ! value) may vary.
456 ! M_A (global) DESCA( M_ ) The number of rows in the global
457 ! array A.
458 ! N_A (global) DESCA( N_ ) The number of columns in the global
459 ! array A.
460 ! MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
461 ! the rows of the array.
462 ! NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
463 ! the columns of the array.
464 ! RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
465 ! row of the array A is distributed.
466 ! CSRC_A (global) DESCA( CSRC_ ) The process column over which the
467 ! first column of the array A is
468 ! distributed.
469 ! LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
470 ! array. LLD_A >= MAX(1,LOCr(M_A)).
471 !
472 ! Let K be the number of rows or columns of a distributed matrix,
473 ! and assume that its process grid has dimension p x q.
474 ! LOCr( K ) denotes the number of elements of K that a process
475 ! would receive if K were distributed over the p processes of its
476 ! process column.
477 ! Similarly, LOCc( K ) denotes the number of elements of K that a
478 ! process would receive if K were distributed over the q processes of
479 ! its process row.
480 ! The values of LOCr() and LOCc() may be determined via a call to the
481 ! ScaLAPACK tool function, NUMROC:
482 ! LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
483 ! LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
484 ! An upper bound for these quantities may be computed by:
485 ! LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
486 ! LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
487 !
488 ! This routine requires square block data decomposition ( MB_A=NB_A ).
489 !
490 ! Arguments
491 ! =========
492 !
493 ! TRANS (global input) CHARACTER
494 ! Specifies the form of the system of equations:
495 ! = 'N': sub( A ) * X = sub( B ) (No transpose)
496 ! = 'T': sub( A )**T * X = sub( B ) (Transpose)
497 ! = 'C': sub( A )**T * X = sub( B ) (Transpose)
498 !
499 ! N (global input) INTEGER
500 ! The number of rows and columns to be operated on, i.e. the
501 ! order of the distributed submatrix sub( A ). N >= 0.
502 !
503 ! NRHS (global input) INTEGER
504 ! The number of right hand sides, i.e., the number of columns
505 ! of the distributed submatrix sub( B ). NRHS >= 0.
506 !
507 ! A (local input) REAL(dp) pointer into the local
508 ! memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
509 ! On entry, this array contains the local pieces of the factors
510 ! L and U from the factorization sub( A ) = P!L!U; the unit
511 ! diagonal elements of L are not stored.
512 !
513 ! IA (global input) INTEGER
514 ! The row index in the global array A indicating the first
515 ! row of sub( A ).
516 !
517 ! JA (global input) INTEGER
518 ! The column index in the global array A indicating the
519 ! first column of sub( A ).
520 !
521 ! DESCA (global and local input) INTEGER array of dimension DLEN_.
522 ! The array descriptor for the distributed matrix A.
523 !
524 ! IPIV (local input) INTEGER array, dimension ( LOCr(M_A)+MB_A )
525 ! This array contains the pivoting information.
526 ! IPIV(i) -> The global row local row i was swapped with.
527 ! This array is tied to the distributed matrix A.
528 !
529 ! B (local input/local output) REAL(dp) pointer into the
530 ! local memory to an array of dimension
531 ! (LLD_B,LOCc(JB+NRHS-1)). On entry, the right hand sides
532 ! sub( B ). On exit, sub( B ) is overwritten by the solution
533 ! distributed matrix X.
534 !
535 ! IB (global input) INTEGER
536 ! The row index in the global array B indicating the first
537 ! row of sub( B ).
538 !
539 ! JB (global input) INTEGER
540 ! The column index in the global array B indicating the
541 ! first column of sub( B ).
542 !
543 ! DESCB (global and local input) INTEGER array of dimension DLEN_.
544 ! The array descriptor for the distributed matrix B.
545 !
546 ! INFO (global output) INTEGER
547 ! = 0: successful exit
548 ! < 0: If the i-th argument is an array and the j-entry had
549 ! an illegal value, then INFO = -(i!100+j), if the i-th
550 ! argument is a scalar and had an illegal value, then
551 ! INFO = -i.
552 !
553 ! =====================================================================
554 !
555 ! .. Parameters ..
556  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, &
557  & LLD_, MB_, M_, NB_, N_, RSRC_
558  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1, &
559  & ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6, &
560  & rsrc_ = 7, csrc_ = 8, lld_ = 9 )
561  REAL(dp) ONE
562  PARAMETER ( ONE = 1.0d+0 )
563 ! ..
564 ! .. Local Scalars ..
565  LOGICAL NOTRAN
566  INTEGER IAROW, IBROW, ICOFFA, ICTXT, IROFFA, IROFFB, &
567  & mycol, myrow, npcol, nprow
568 ! ..
569 ! .. Local Arrays ..
570  INTEGER DESCIP( DLEN_ ), IDUM1( 1 ), IDUM2( 1 )
571 ! ..
572 ! .. External Subroutines ..
573  EXTERNAL blacs_gridinfo, chk1mat, descset, pchk2mat, &
574  & pdlapiv, pdtrsm, pxerbla
575 ! ..
576 ! .. External Functions ..
577  LOGICAL LSAME
578  INTEGER INDXG2P, NUMROC
579  EXTERNAL indxg2p, lsame, numroc
580 ! ..
581 ! .. Intrinsic Functions ..
582  INTRINSIC ichar, mod
583 ! ..
584 ! .. Executable Statements ..
585 !
586 ! Get grid parameters
587 !
588  ictxt = desca( ctxt_ )
589  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
590 !
591 ! Test the input parameters
592 !
593  info = 0
594  IF( nprow.EQ.-1 ) THEN
595  info = -(700+ctxt_)
596  ELSE
597  notran = lsame( trans, 'N' )
598  CALL chk1mat( n, 2, n, 2, ia, ja, desca, 7, info )
599  CALL chk1mat( n, 2, nrhs, 3, ib, jb, descb, 12, info )
600  IF( info.EQ.0 ) THEN
601  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ), &
602  & nprow )
603  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ), &
604  & nprow )
605  iroffa = mod( ia-1, desca( mb_ ) )
606  icoffa = mod( ja-1, desca( nb_ ) )
607  iroffb = mod( ib-1, descb( mb_ ) )
608  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT. &
609  & lsame( trans, 'C' ) ) THEN
610  info = -1
611  ELSE IF( iroffa.NE.0 ) THEN
612  info = -5
613  ELSE IF( icoffa.NE.0 ) THEN
614  info = -6
615  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
616  info = -(700+nb_)
617 ! ELSE IF( IROFFB.NE.0 .OR. IBROW.NE.IAROW ) THEN
618 ! INFO = -10
619  ELSE IF( descb( mb_ ).NE.desca( nb_ ) ) THEN
620  info = -(1200+nb_)
621  ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
622  info = -(1200+ctxt_)
623  END IF
624  END IF
625  IF( notran ) THEN
626  idum1( 1 ) = ichar( 'N' )
627  ELSE IF( lsame( trans, 'T' ) ) THEN
628  idum1( 1 ) = ichar( 'T' )
629  ELSE
630  idum1( 1 ) = ichar( 'C' )
631  END IF
632  idum2( 1 ) = 1
633  CALL pchk2mat( n, 2, n, 2, ia, ja, desca, 7, n, 2, nrhs, 3, &
634  & ib, jb, descb, 12, 1, idum1, idum2, info )
635  END IF
636 
637  IF( info.NE.0 ) THEN
638  CALL pxerbla( ictxt, 'PDGETRS', -info )
639  RETURN
640  END IF
641 !
642 ! Quick return if possible
643 !
644  IF( n.EQ.0 .OR. nrhs.EQ.0 ) RETURN
645 
646  CALL descset( descip, desca( m_ ) + desca( mb_ )*nprow, 1, &
647  & desca( mb_ ), 1, desca( rsrc_ ), mycol, ictxt, &
648  & desca( mb_ ) + numroc( desca( m_ ), desca( mb_ ), &
649  & myrow, desca( rsrc_ ), nprow ) )
650 
651  IF( notran ) THEN
652 !
653 ! Solve sub( A ) * X = sub( B ).
654 !
655 ! Apply row interchanges to the right hand sides.
656 !
657  CALL pdlapiv( 'Forward', 'Row', 'Col', n, nrhs, b, ib, jb, &
658  & descb, ipiv, ia, 1, descip, idum1 )
659 !
660 ! Solve L*X = sub( B ), overwriting sub( B ) with X.
661 !
662  CALL pdtrsm( 'Left', 'Lower', 'No transpose', 'Unit', n, nrhs, &
663  & one, a, ia, ja, desca, b, ib, jb, descb )
664 !
665 ! Solve U*X = sub( B ), overwriting sub( B ) with X.
666 !
667  CALL pdtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', n, &
668  & nrhs, one, a, ia, ja, desca, b, ib, jb, descb )
669  ELSE
670 !
671 ! Solve sub( A )' * X = sub( B ).
672 !
673 ! Solve U'*X = sub( B ), overwriting sub( B ) with X.
674 !
675  CALL pdtrsm( 'Left', 'Upper', 'Transpose', 'Non-unit', n, nrhs, &
676  & one, a, ia, ja, desca, b, ib, jb, descb )
677 !
678 ! Solve L'*X = sub( B ), overwriting sub( B ) with X.
679 !
680  CALL pdtrsm( 'Left', 'Lower', 'Transpose', 'Unit', n, nrhs, &
681  & one, a, ia, ja, desca, b, ib, jb, descb )
682 !
683 ! Apply row interchanges to the solution vectors.
684 !
685  CALL pdlapiv( 'Backward', 'Row', 'Col', n, nrhs, b, ib, jb, &
686  & descb, ipiv, ia, 1, descip, idum1 )
687 
688  END IF
689 
690  RETURN
691 
692 ! End of PDGETRS
693 
694  END SUBROUTINE pdgetrs
695 
696  SUBROUTINE pdelget( SCOPE, TOP, ALPHA, A, IA, JA, DESCA )
697  use descriptor_mod, desca_x=>desca
698  implicit none
699 !
700 ! -- ScaLAPACK tools routine (version 1.7) --
701 ! University of Tennessee, Knoxville, Oak Ridge National Laboratory,
702 ! and University of California, Berkeley.
703 ! May 1, 1997
704 !
705 ! .. Scalar Arguments ..
706  CHARACTER*1 SCOPE, TOP
707  INTEGER IA, JA
708  REAL(dp) ALPHA
709 ! ..
710 ! .. Array arguments ..
711  INTEGER DESCA( * )
712  REAL(dp) A( * )
713 ! ..
714 !
715 ! Purpose
716 ! =======
717 !
718 ! PDELGET sets alpha to the distributed matrix entry A( IA, JA ).
719 ! The value of alpha is set according to the scope.
720 !
721 ! Notes
722 ! =====
723 !
724 ! Each global data object is described by an associated description
725 ! vector. This vector stores the information required to establish
726 ! the mapping between an object element and its corresponding process
727 ! and memory location.
728 !
729 ! Let A be a generic term for any 2D block cyclicly distributed array.
730 ! Such a global array has an associated description vector DESCA.
731 ! In the following comments, the character _ should be read as
732 ! "of the global array".
733 !
734 ! NOTATION STORED IN EXPLANATION
735 ! --------------- -------------- --------------------------------------
736 ! DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
737 ! DTYPE_A = 1.
738 ! CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
739 ! the BLACS process grid A is distribu-
740 ! ted over. The context itself is glo-
741 ! bal, but the handle (the integer
742 ! value) may vary.
743 ! M_A (global) DESCA( M_ ) The number of rows in the global
744 ! array A.
745 ! N_A (global) DESCA( N_ ) The number of columns in the global
746 ! array A.
747 ! MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
748 ! the rows of the array.
749 ! NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
750 ! the columns of the array.
751 ! RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
752 ! row of the array A is distributed.
753 ! CSRC_A (global) DESCA( CSRC_ ) The process column over which the
754 ! first column of the array A is
755 ! distributed.
756 ! LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
757 ! array. LLD_A >= MAX(1,LOCr(M_A)).
758 !
759 ! Let K be the number of rows or columns of a distributed matrix,
760 ! and assume that its process grid has dimension p x q.
761 ! LOCr( K ) denotes the number of elements of K that a process
762 ! would receive if K were distributed over the p processes of its
763 ! process column.
764 ! Similarly, LOCc( K ) denotes the number of elements of K that a
765 ! process would receive if K were distributed over the q processes of
766 ! its process row.
767 ! The values of LOCr() and LOCc() may be determined via a call to the
768 ! ScaLAPACK tool function, NUMROC:
769 ! LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
770 ! LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
771 ! An upper bound for these quantities may be computed by:
772 ! LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )!MB_A
773 ! LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )!NB_A
774 !
775 ! Arguments
776 ! =========
777 !
778 ! SCOPE (global input) CHARACTER*1
779 ! The BLACS scope in which alpha is updated.
780 ! If SCOPE = 'R', alpha is updated only in the process row
781 ! containing A( IA, JA ),
782 ! If SCOPE = 'C', alpha is updated only in the process column
783 ! containing A( IA, JA ),
784 ! If SCOPE = 'A', alpha is updated in all the processes of the
785 ! grid,
786 ! otherwise alpha is updated only in the process containing
787 ! A( IA, JA ).
788 !
789 ! TOP (global input) CHARACTER!1
790 ! The topology to be used if broadcast is needed.
791 !
792 ! ALPHA (global output) DOUBLE PRECISION, the scalar alpha.
793 !
794 ! A (local input) REAL(dp) pointer into the local memory
795 ! to an array of dimension (LLD_A,!) containing the local
796 ! pieces of the distributed matrix A.
797 !
798 ! IA (global input) INTEGER
799 ! The row index in the global array A indicating the first
800 ! row of sub( A ).
801 !
802 ! JA (global input) INTEGER
803 ! The column index in the global array A indicating the
804 ! first column of sub( A ).
805 !
806 ! DESCA (global and local input) INTEGER array of dimension DLEN_.
807 ! The array descriptor for the distributed matrix A.
808 !
809 ! =====================================================================
810 !
811 ! .. Parameters ..
812 ! INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
813 ! $ LLD_, MB_, M_, NB_, N_, RSRC_
814 ! PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
815 ! $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
816 ! $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
817  REAL(dp) ZERO
818  PARAMETER ( ZERO = 0.0d+0 )
819 ! ..
820 ! .. Local Scalars ..
821  INTEGER IACOL, IAROW, ICTXT, IIA, IOFFA, JJA !, MYCOL, MYROW, NPCOL, NPROW
822 ! ..
823 ! .. External Subroutines ..
824  EXTERNAL blacs_gridinfo, dgebr2d, dgebs2d, infog2l
825 ! ..
826 ! .. External Functions ..
827  LOGICAL LSAME
828  EXTERNAL LSAME
829 ! ..
830 ! .. Executable Statements ..
831 !
832 ! Get grid parameters.
833 !
834  ictxt = desca( ctxt_ )
835  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
836 !
837  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,&
838  iarow, iacol )
839 !
840  alpha = zero
841 !
842  IF( lsame( scope, 'R' ) ) THEN
843  IF( myrow.EQ.iarow ) THEN
844  IF( mycol.EQ.iacol ) THEN
845  ioffa = iia+(jja-1)*desca( lld_ )
846  CALL dgebs2d( ictxt, scope, top, 1, 1, a( ioffa ), 1 )
847  alpha = a( ioffa )
848  ELSE
849  CALL dgebr2d( ictxt, scope, top, 1, 1, alpha, 1, &
850  iarow, iacol )
851  END IF
852  END IF
853  ELSE IF( lsame( scope, 'C' ) ) THEN
854  IF( mycol.EQ.iacol ) THEN
855  IF( myrow.EQ.iarow ) THEN
856  ioffa = iia+(jja-1)*desca( lld_ )
857  CALL dgebs2d( ictxt, scope, top, 1, 1, a( ioffa ), 1 )
858  alpha = a( ioffa )
859  ELSE
860  CALL dgebr2d( ictxt, scope, top, 1, 1, alpha, 1, &
861  iarow, iacol )
862  END IF
863  END IF
864  ELSE IF( lsame( scope, 'A' ) ) THEN
865  IF( ( myrow.EQ.iarow ).AND.( mycol.EQ.iacol ) ) THEN
866  ioffa = iia+(jja-1)*desca( lld_ )
867  CALL dgebs2d( ictxt, scope, top, 1, 1, a( ioffa ), 1 )
868  alpha = a( ioffa )
869  ELSE
870  CALL dgebr2d( ictxt, scope, top, 1, 1, alpha, 1, &
871  iarow, iacol )
872  END IF
873  ELSE
874  IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) &
875  alpha = a( iia+(jja-1)*desca( lld_ ) )
876  END IF
877 !
878  RETURN
879 !
880 ! End of PDELGET
881 !
882  END
883 
884  SUBROUTINE pdelset( A, IA, JA, DESCA, ALPHA )
885  use descriptor_mod, desca_x=>desca
886  implicit none
887 !
888 ! -- ScaLAPACK tools routine (version 1.7) --
889 ! University of Tennessee, Knoxville, Oak Ridge National Laboratory,
890 ! and University of California, Berkeley.
891 ! May 1, 1997
892 !
893 ! .. Scalar Arguments ..
894  INTEGER IA, JA
895  REAL(dp) ALPHA
896 ! ..
897 ! .. Array arguments ..
898  INTEGER DESCA( * )
899  REAL(dp) A( * )
900 ! ..
901 !
902 ! Purpose
903 ! =======
904 !
905 ! PDELSET sets the distributed matrix entry A( IA, JA ) to ALPHA.
906 !
907 ! Notes
908 ! =====
909 !
910 ! Each global data object is described by an associated description
911 ! vector. This vector stores the information required to establish
912 ! the mapping between an object element and its corresponding process
913 ! and memory location.
914 !
915 ! Let A be a generic term for any 2D block cyclicly distributed array.
916 ! Such a global array has an associated description vector DESCA.
917 ! In the following comments, the character _ should be read as
918 ! "of the global array".
919 !
920 ! NOTATION STORED IN EXPLANATION
921 ! --------------- -------------- --------------------------------------
922 ! DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
923 ! DTYPE_A = 1.
924 ! CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
925 ! the BLACS process grid A is distribu-
926 ! ted over. The context itself is glo-
927 ! bal, but the handle (the integer
928 ! value) may vary.
929 ! M_A (global) DESCA( M_ ) The number of rows in the global
930 ! array A.
931 ! N_A (global) DESCA( N_ ) The number of columns in the global
932 ! array A.
933 ! MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
934 ! the rows of the array.
935 ! NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
936 ! the columns of the array.
937 ! RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
938 ! row of the array A is distributed.
939 ! CSRC_A (global) DESCA( CSRC_ ) The process column over which the
940 ! first column of the array A is
941 ! distributed.
942 ! LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
943 ! array. LLD_A >= MAX(1,LOCr(M_A)).
944 !
945 ! Let K be the number of rows or columns of a distributed matrix,
946 ! and assume that its process grid has dimension p x q.
947 ! LOCr( K ) denotes the number of elements of K that a process
948 ! would receive if K were distributed over the p processes of its
949 ! process column.
950 ! Similarly, LOCc( K ) denotes the number of elements of K that a
951 ! process would receive if K were distributed over the q processes of
952 ! its process row.
953 ! The values of LOCr() and LOCc() may be determined via a call to the
954 ! ScaLAPACK tool function, NUMROC:
955 ! LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
956 ! LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
957 ! An upper bound for these quantities may be computed by:
958 ! LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
959 ! LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
960 !
961 ! Arguments
962 ! =========
963 !
964 ! A (local output) REAL(dp) pointer into the local memory
965 ! to an array of dimension (LLD_A,*) containing the local
966 ! pieces of the distributed matrix A.
967 !
968 ! IA (global input) INTEGER
969 ! The row index in the global array A indicating the first
970 ! row of sub( A ).
971 !
972 ! JA (global input) INTEGER
973 ! The column index in the global array A indicating the
974 ! first column of sub( A ).
975 !
976 ! DESCA (global and local input) INTEGER array of dimension DLEN_.
977 ! The array descriptor for the distributed matrix A.
978 !
979 ! ALPHA (local input) DOUBLE PRECISION
980 ! The scalar alpha.
981 !
982 ! =====================================================================
983 !
984 ! .. Parameters ..
985 ! INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
986 ! $ LLD_, MB_, M_, NB_, N_, RSRC_
987 ! PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
988 ! $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
989 ! $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
990 ! ..
991 ! .. Local Scalars ..
992  INTEGER IACOL, IAROW, IIA, JJA !, MYCOL, MYROW, NPCOL, NPROW
993 ! ..
994 ! .. External Subroutines ..
995  EXTERNAL blacs_gridinfo, infog2l
996 ! ..
997 ! .. Executable Statements ..
998 !
999 ! Get grid parameters.
1000 !
1001  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
1002 !
1003  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,&
1004  iarow, iacol )
1005 !
1006  IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) &
1007  a( iia+(jja-1)*desca( lld_ ) ) = alpha
1008 !
1009  RETURN
1010 !
1011 ! End of PDELSET
1012 !
1013  END
1014 #endif
1015 
1016 #endif
1017 
1018  end module ptrd_mod
v3_utilities::assert
Definition: v3_utilities.f:55