New swap data code.

psblas3-type-indexed
Salvatore Filippone 20 years ago
parent 2d87bbae3b
commit 30dc67fbc0

@ -84,7 +84,7 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
ix = 1 ix = 1
if (present(jx)) then if (present(jx)) then
ijx = jx ijx = jx
@ -184,7 +184,7 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
end if end if
if(info.ne.0) then if(info.ne.0) then
ch_err='PSI_dSwap...' ch_err='PSI_dSwapdata'
call psb_errpush(4010,name,a_err=ch_err) call psb_errpush(4010,name,a_err=ch_err)
goto 9999 goto 9999
end if end if

@ -1,344 +0,0 @@
/*
* Parallel Sparse BLAS v2.0
* (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
* Alfredo Buttari University of Rome Tor Vergata
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions, and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. The name of the PSBLAS group or the names of its contributors may
* not be used to endorse or promote products derived from this
* software without specific written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*
*/
/* This header file replaces every call to a BLACS routine by C interface
with the same call performed by Fortran interface */
#ifndef CTOF_BLACS
#define CTOF_BLACS
#endif
/* Variables necessary for invocations where
constant arguments are used */
static int i1, i2, i3, i4, i5, i6, i7;
/* Support routines:
Initialization */
#define Cblacs_pinfo(mypnum, nprocs) \
blacs_pinfo_(mypnum, nprocs)
#define Cblacs_setup(mypnum, nprocs) \
blacs_setup_(mypnum, nprocs)
#define Cblacs_get(ictxt, what, val) \
{i1 = ictxt; i2 = what; \
blacs_get_(&i1, &i2,val);}
#define Cblacs_set(ictxt, what, val) \
{i1 = ictxt; i2 = what; \
blacs_set_(&i1, &i2, &val);}
#define Cblacs_gridinit(ictxt, order, nprow, npcol) \
{i1 = nprow; i2 = npcol; \
blacs_gridinit_(ictxt, order, &i1, &i2);}
#define Cblacs_gridmap(ictxt, pmap, ldpmap, nprow, npcol) \
{i1 = ldpmap; i2 = nprow; i3 = npcol; \
blacs_gridmap_(ictxt, pmap, &i1, &i2, &i3);}
/* Support routines:
Destruction */
#define Cblacs_freebuff(ictxt, wait) \
{i1 = ictxt; i2 = wait; \
blacs_freebuff_(&i1, &i2);}
#define Cblacs_gridexit(ictxt) \
{i1 = ictxt; \
blacs_gridexit_(&i1);}
#define Cblacs_abort(ictxt, errornum) \
{i1 = ictxt; i2 = errornum; \
blacs_abort_(&i1, &i2);}
#define Cblacs_exit(doneflag) \
{i1 = doneflag; \
blacs_exit_(&i1);}
/* Support routines:
Informational and Miscellaneous */
#define Cblacs_gridinfo(ictxt,nprow,npcol,myprow,mypcol) \
{i1 = ictxt; \
blacs_gridinfo_(&i1, nprow, npcol, myprow, mypcol);}
#define Cblacs_pnum(ictxt, prow, pcol) \
{i1 = ictxt; i2 = prow; i3 = pcol; \
blacs_pnum_(&i1, &i2, &i3);}
#define Cblacs_pcoord(ictxt, pnum, prow, pcol) \
{i1 = ictxt; i2 = pnum; \
blacs_pcoord_(&i1, &i2, prow, pcol);}
#define Cblacs_barrier(ictxt, scope) \
{i1 = ictxt; \
blacs_barrier_(&i1, scope);}
/* Support routines:
Unofficial */
#define Csetpvmtids(ntasks, tids) \
{i1 = ntasks; \
setpvmtids_(&i1, tids);}
#define Cdcputime() \
dcputime_()
#define Cdwalltime() \
dwalltime_()
#define Cksendid(ictxt, rdest, cdest) \
{i1 = ictxt; i2 = rdest; i3 = cdest; \
ksendid_(&i1, &i2, &i3);}
#define Ckrecvid(ictxt, rsrc, csrc) \
{i1 = ictxt; i2 = rsrc; i3 = csrc; \
krecvid_(&i1, &i2, &i3);}
#define Ckbsid(ictxt, scope) \
{i1 = ictxt; \
kbsid_(&i1, scope);}
#define Ckbrid(ictxt, scope, rsrc, csrc) \
{i1 = ictxt; i2 = rsrc; i3 = csrc; \
kbrid_(&i1, scope, &i2, &i3);}
/* Point to Point :
Integer */
#define Cigesd2d(ictxt, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
igesd2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Cigerv2d(ictxt, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
igerv2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Citrsd2d(ictxt, uplo, diag, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
itrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
#define Citrrv2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
itrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Point to Point :
Single precision real */
#define Csgesd2d(ictxt, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
sgesd2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Csgerv2d(ictxt, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
sgerv2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Cstrsd2d(ictxt, uplo, diag, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
strsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
#define Cstrrv2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
strsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Point to Point :
Double precision real */
#define Cdgesd2d(ictxt, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
dgesd2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Cdgerv2d(ictxt, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
dgerv2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Cdtrsd2d(ictxt, uplo, diag, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
dtrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
#define Cdtrrv2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
dtrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Point to Point :
Single precision complex */
#define Ccgesd2d(ictxt, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
cgesd2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Ccgerv2d(ictxt, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
cgerv2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Cctrsd2d(ictxt, uplo, diag, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
ctrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
#define Cctrrv2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
ctrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Point to Point :
Double precision complex */
#define Czgesd2d(ictxt, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
zgesd2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Czgerv2d(ictxt, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
zgerv2d_(&i1, &i2, &i3, A, &i4, &i5, &i6);}
#define Cztrsd2d(ictxt, uplo, diag, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
ztrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
#define Cztrrv2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
ztrsd2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Broadcasts :
Integer */
#define Cigebs2d(ictxt, scope, top, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
igebs2d_(&i1, scope, top, &i2, &i3, A, &i4);}
#define Cigebr2d(ictxt, scope, top, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
igebr2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Citrbs2d(ictxt, scope, top, uplo, diag, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
itrbs2d_(&i1, scope, top, uplo, diag, &i2, &i3, A, &i4);}
#define Citrbr2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
igebr2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Broadcasts :
Single precision real */
#define Csgebs2d(ictxt, scope, top, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
sgebs2d_(&i1, scope, top, &i2, &i3, A, &i4);}
#define Csgebr2d(ictxt, scope, top, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
sgebr2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Cstrbs2d(ictxt, scope, top, uplo, diag, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
strbs2d_(&i1, scope, top, uplo, diag, &i2, &i3, A, &i4);}
#define Cstrbr2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
sgebr2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Broadcasts :
Double precision real */
#define Cdgebs2d(ictxt, scope, top, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
dgebs2d_(&i1, scope, top, &i2, &i3, A, &i4);}
#define Cdgebr2d(ictxt, scope, top, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
dgebr2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Cdtrbs2d(ictxt, scope, top, uplo, diag, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
dtrbs2d_(&i1, scope, top, uplo, diag, &i2, &i3, A, &i4);}
#define Cdtrbr2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
dgebr2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Broadcasts :
Single precision complex */
#define Ccgebs2d(ictxt, scope, top, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
cgebs2d_(&i1, scope, top, &i2, &i3, A, &i4);}
#define Ccgebr2d(ictxt, scope, top, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
cgebr2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Cctrbs2d(ictxt, scope, top, uplo, diag, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
ctrbs2d_(&i1, scope, top, uplo, diag, &i2, &i3, A, &i4);}
#define Cctrbr2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
cgebr2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Broadcasts :
Double precision complex */
#define Czgebs2d(ictxt, scope, top, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
zgebs2d_(&i1, scope, top, &i2, &i3, A, &i4);}
#define Czgebr2d(ictxt, scope, top, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
zgebr2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Cztrbs2d(ictxt, scope, top, uplo, diag, m, n, A, lda) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; \
ztrbs2d_(&i1, scope, top, uplo, diag, &i2, &i3, A, &i4);}
#define Cztrbr2d(ictxt, uplo, diag, m, n, A, lda, rsrc, csrc) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rsrc; i6 = csrc; \
zgebr2d_(&i1, uplo, diag, &i2, &i3, A, &i4, &i5, &i6);}
/* Combines:
Integer */
#define Cigsum2d(ictxt, scope, top, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
igsum2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Cigamx2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
igamx2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
#define Cigamn2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
igamn2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
/* Combines:
Single precision real */
#define Csgsum2d(ictxt, scope, top, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
sgsum2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Csgamx2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
sgamx2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
#define Csgamn2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
sgamn2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
/* Combines:
Double precision real */
#define Cdgsum2d(ictxt, scope, top, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
dgsum2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Cdgamx2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
dgamx2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
#define Cdgamn2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
dgamn2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
/* Combines:
Single precision complex */
#define Ccgsum2d(ictxt, scope, top, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
cgsum2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Ccgamx2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
cgamx2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
#define Ccgamn2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
cgamn2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
/* Combines:
Double precision complex */
#define Czgsum2d(ictxt, scope, top, m, n, A, lda, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = rdest; i6 = cdest; \
zgsum2d_(&i1, scope, top, &i2, &i3, A, &i4, &i5, &i6);}
#define Czgamx2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
zgamx2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}
#define Czgamn2d(ictxt, scope, top, m, n, A, lda, RA, CA, RCflag, rdest, cdest) \
{i1 = ictxt; i2 = m; i3 = n; i4 = lda; i5 = RCflag; i6 = rdest; i7 = cdest; \
zgamn2d_(&i1, scope, top, &i2, &i3, A, &i4, RA, CA, &i5, &i6, &i7);}

@ -1,316 +0,0 @@
/*
* Parallel Sparse BLAS v2.0
* (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
* Alfredo Buttari University of Rome Tor Vergata
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions, and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. The name of the PSBLAS group or the names of its contributors may
* not be used to endorse or promote products derived from this
* software without specific written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*
*/
/*
* This file includes the standard C libraries, as well as system
* dependent include files. All PSBLAS routines include this file.
*/
#include <string.h>
#ifndef PSBLASH
#define PSBLASH
/*
* ========================================================================
* Machine Specific PBLAS macros
* ========================================================================
*/
/* This is a debugging option.
#define PS_CONTROL_LEVEL */
#define _HAL_ 0
#define _T3D_ 1
#ifdef T3D
#define _MACH_ _T3D_
#endif
#ifndef _MACH_
#define _MACH_ _HAL_
#endif
/*
* ========================================================================
* Include files
* ========================================================================
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#if( _MACH_ == _T3D_ )
#include <fortran.h>
#endif
#ifdef USE_FBLACS
#ifndef CTOF_BLACS
#include "ctof_blacs.h"
#endif
#endif
/*
* ========================================================================
* FORTRAN <-> C interface
* ========================================================================
*
* These macros define how the PBLAS will be called. _F2C_ADD_ assumes
* that they will be called by FORTRAN, which expects C routines to have
* an underscore postfixed to the name (Suns, and Intel machines expect
* this). _F2C_NOCHANGE indicates that FORTRAN will be calling, and that
* it expects the name called by FORTRAN to be identical to that compiled
* by the C (RS6K's do this). _F2C_UPCASE says it expects C routines
* called by FORTRAN to be in all upcase (CRAY wants this).
*/
#define _F2C_ADD_ 0
#define _F2C_NOCHANGE 1
#define _F2C_UPCASE 2
#ifdef UpCase
#define _F2C_CALL_ _F2C_UPCASE
#endif
#ifdef NoChange
#define _F2C_CALL_ _F2C_NOCHANGE
#endif
#ifdef Add_
#define _F2C_CALL_ _F2C_ADD_
#endif
#ifndef _F2C_CALL_
#define _F2C_CALL_ _F2C_ADD_
#endif
/*
* ========================================================================
* TYPE DEFINITIONS AND CONVERSION UTILITIES
* ========================================================================
*/
typedef struct { float re, im; } complex;
typedef struct { double re, im; } complex16;
#if( _MACH_ == _T3D_ )
/* Type of character argument in a FORTRAN call */
#define F_CHAR _fcd
/* Character conversion utilities */
#define F2C_CHAR(a) ( _fcdtocp( (a) ) )
#define C2F_CHAR(a) ( _cptofcd( (a), 1 ) )
/* Type of FORTRAN functions */
#define F_VOID_FCT void fortran /* Subroutine */
#define F_INTG_FCT int fortran /* INTEGER function */
#define F_DBLE_FCT double fortran /* DOUBLE PRECISION function */
#else
/* Type of character argument in a FORTRAN call */
typedef char * F_CHAR;
/* Character conversion utilities */
#define F2C_CHAR(a) (a)
#define C2F_CHAR(a) (a)
/* Type of FORTRAN functions */
#define F_VOID_FCT void /* Subroutine */
#define F_INTG_FCT int /* INTEGER function */
#define F_DBLE_FCT double /* DOUBLE PRECISION function */
#endif
/*
* ======================================================================
* FUNCTIONS PROTOTYPES
* ======================================================================
*/
void DVSct(int n, int k,int idx[],int flag, double X[], int lx,
double beta, double Y[], int ly);
void DVGth(int n, int k,int idx[],int flag, double X[], int lx,double Y[], int ly);
void IVSct(int n, int k,int idx[],int flag, int X[], int lx,
int beta, int Y[], int ly);
void IVGth(int n, int k,int idx[],int flag, int X[], int lx,int Y[], int ly);
void PSI_dSwapData(int iflag, int n, double beta, double Y[], int ly,
int desc_data[], int desc_halo[],
double *work, int *lwork, int *ierror);
void PSI_dSwapTran(int flag, int n, double beta, double Y[], int ly,
int desc_data[], int desc_halo[],
double *work, int *lwork, int *ierror);
void PSI_zSwapData(int n, double Y[], int ly, int desc_data[], int desc_halo[],
double *work, int *lwork, int *ierror);
void PSI_zSwapOverlap(double Y[], double Sum_Ovrlap[], int desc_data[],
int desc_ovrlap[], double work[], int *lwork, int *ierror);
void PSI_iSwapData(int iflag, int n, int beta, int Y[], int ly,
int desc_data[], int desc_halo[],
int *work, int *lwork, int *ierror);
void PSI_iSwapTran(int flag, int n, int beta, int Y[], int ly,
int desc_data[], int desc_halo[],
int *work, int *lwork, int *ierror);
/*
* ========================================================================
* #DEFINE MACRO CONSTANTS
* ========================================================================
*/
/* MACRO max */
#define max(x,y) ((x)>(y)?(x):(y))
/*MACRO for ovrlap update*/
#define NOHALO_ 0
#define HALO_ 4
#define NONE_ 0
#define SUM_ 1
#define AVG_ 2
#define SQUARE_ROOT_ 3
/* Bit fields to control swapdata/ovrlap behaviour.
BEWARE: check consistency with tools_const.f.
Should it be automated? */
#define SWAP_SEND 1
#define SWAP_RECV 2
#define SWAP_SYNC 4
#define SWAP_MPI 8
/* Macro for MATRIX_DATA array */
#define DEC_TYPE_ 0 /* The type of decomposition of global
matrix A. */
#define M_ 1 /* Number of equations */
#define N_ 2 /* Number of variables */
#define N_ROW_ 3 /* The number of row of local matrix. */
#define N_COL_ 4 /* The number of columns of local
matrix. */
#define CTXT_ 5 /* The BLACS context handle, indicating
the global context of the operation
on the matrix.
The context itself is global. */
#define LOC_TO_GLOB_ 6 /* The pointer to the array
loc_to_glob */
#define MPI_C_ 8 /* The MPI Fortran handle */
/* values for DEC_TYPE_ */
#define DESC_ASB 3099
#define DESC_BLD (DESC_ASB+1)
/* Macro for HALO array */
#define PROC_ID_ 0 /* The identifier of domain. */
#define N_ELEM_RECV_ 1 /* The number of elements to receive*/
#define ELEM_RECV_ 2 /* The first index of local elements */
#define N_ELEM_SEND_ 2 /* The number of elements to send */
#define ELEM_SEND_ 3 /* The first index of local elements */
/* Macro for OVERLAP array */
#define N_OVRLP_ELEM_ 1 /* The number of overlap elements to recv/send */
#define OVRLP_ELEM_TO_ 2 /* The first index of local elements */
/* Macro for OVR_ELEM_D array */
#define OVRLP_ELEM_ 0
#define N_DOM_OVR_ 1
#define BROADCAST "B" /* Blacs operation definitions */
#define COMBINE "C"
#define ALL "A" /* Scope definitions */
#define COLUMN "C"
#define ROW "R"
#define TOPDEF " " /* Default BLACS topology, PB-BLAS routines */
#define CTOPDEF ' '
#define TOPGET "!"
#define YES "Y"
#define NO "N"
#define MULLENFAC 2
#define ONE 1.0
#define ZERO 0.0
/* Integer values for error checking */
#define no_err 0
#define act_ret 0
#define act_abort 1
/*
* ========================================================================
* PREPROCESSOR MACRO FUNCTIONS USED FOR OPTIMIZATION & CONVENIENCE
* ========================================================================
*/
#define ABS(a) ((a > 0) ? (a) : (-a))
#define MIN(a,b) ((a < b) ? (a) : (b))
#define MAX(a,b) ((a > b) ? (a) : (b))
#define CEIL(a,b) ( (a+b-1) / (b) )
#define Mlowcase(C) ( ((C) > 64 && (C) < 91) ? (C) | 32 : (C) )
#define Mupcase(C) ( ((C) > 96 && (C) < 123) ? (C) & 0xDF : (C) )
#define INDXG2L( iglob, nb, iproc, isrcproc, nprocs )\
( (nb) * ( ( (iglob)-1) / ( (nb) * (nprocs) ) ) +\
( ( (iglob) - 1 ) % (nb) ) + 1 )
#define INDXL2G( iloc, nb, iproc, isrcproc, nprocs )\
( (nprocs) * (nb) * ( ( (iloc) - 1 ) / (nb) ) +\
( ( (iloc) - 1 ) % (nb) ) +\
( ( (nprocs) + (iproc) - (isrcproc) ) % (nprocs) ) * (nb) + 1 )
#define INDXG2P( iglob, nb, iproc, isrcproc, nprocs ) \
( ( (isrcproc) + ( (iglob) - 1 ) / (nb) ) % (nprocs) )
#define MYROC0( nblocks, n, nb, nprocs )\
( ( (nblocks) % (nprocs) ) ? ( ( (nblocks) / (nprocs) ) * (nb) + (nb) )\
: ( ( (nblocks) / (nprocs) )* (nb) + ( (n) % (nb) ) ) )
#if( _F2C_CALL_ == _F2C_ADD_ )
/*
* These defines set up the naming scheme required to have a FORTRAN
* routine call a C routine (which is what the PBLAS are written in).
* No redefinition necessary to have following FORTRAN to C interface:
* FORTRAN CALL C DECLARATION
* call pdgemm(...) void pdgemm_(...)
*
* This is the default.
*/
#define pbchkvectf pbchkvectf_
#define fcpsb_errcomm fcpsb_errcomm_
#define fcpsb_erractionsave fcpsb_erractionsave_
#define fcpsb_erractionrestore fcpsb_erractionrestore_
#define fcpsb_perror fcpsb_perror_
#define fcpsb_serror fcpsb_serror_
#define fcpsb_errpush fcpsb_errpush_
#endif

@ -28,7 +28,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info) subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info)
use psb_realloc_mod use psb_realloc_mod
use psb_descriptor_type use psb_descriptor_type
@ -37,7 +37,7 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info)
implicit none implicit none
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info integer, intent(out) :: info,nxch,nsnd,nrcv
integer, intent(in) :: index_in(:) integer, intent(in) :: index_in(:)
integer, pointer :: index_out(:) integer, pointer :: index_out(:)
logical :: glob_idx logical :: glob_idx
@ -50,13 +50,6 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info)
logical,parameter :: debug=.false. logical,parameter :: debug=.false.
character(len=20) :: name character(len=20) :: name
interface
subroutine psi_compute_size(desc_data,&
& index_in, dl_lda, info)
integer :: info, dl_lda
integer :: desc_data(:), index_in(:)
end subroutine psi_compute_size
end interface
interface interface
subroutine psi_sort_dl(dep_list,l_dep_list,np,info) subroutine psi_sort_dl(dep_list,l_dep_list,np,info)
@ -73,12 +66,12 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info)
interface interface
subroutine psi_desc_index(desc_data,index_in,dep_list,& subroutine psi_desc_index(desc_data,index_in,dep_list,&
& length_dl,loc_to_glob,glob_to_loc,desc_index,& & length_dl,nsnd,nrcv,loc_to_glob,glob_to_loc,desc_index,&
& isglob_in,info) & isglob_in,info)
integer :: desc_data(:),index_in(:),dep_list(:) integer :: desc_data(:),index_in(:),dep_list(:)
integer :: loc_to_glob(:),glob_to_loc(:) integer :: loc_to_glob(:),glob_to_loc(:)
integer,pointer :: desc_index(:) integer,pointer :: desc_index(:)
integer :: length_dl, info integer :: length_dl,nsnd,nrcv, info
logical :: isglob_in logical :: isglob_in
end subroutine psi_desc_index end subroutine psi_desc_index
end interface end interface
@ -96,7 +89,6 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info)
endif endif
! allocate dependency list ! allocate dependency list
! call psi_compute_size(desc_a%matrix_data, index_in, dl_lda, info)
! This should be computed more efficiently to save space when ! This should be computed more efficiently to save space when
! the number of processors becomes very high ! the number of processors becomes very high
dl_lda=np+1 dl_lda=np+1
@ -119,10 +111,10 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info)
goto 9999 goto 9999
end if end if
if (debug) write(*,*) 'crea_index: from extract_dep_list',& if (debug) write(0,*) 'crea_index: from extract_dep_list',&
& me,length_dl(0),index_in(1), ':',dep_list(:length_dl(me),me) & me,length_dl(0),index_in(1), ':',dep_list(:length_dl(me),me)
! ...now process root contains dependence list of all processes... ! ...now process root contains dependence list of all processes...
if (debug) write(*,*) 'crea_halo: root sorting dep list' if (debug) write(0,*) 'crea_index: root sorting dep list'
! ....i must order communication in in halo ! ....i must order communication in in halo
call psi_dl_check(dep_list,max(1,dl_lda),np,length_dl) call psi_dl_check(dep_list,max(1,dl_lda),np,length_dl)
@ -138,9 +130,11 @@ subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info)
if(debug) write(0,*)'in psi_crea_index calling psi_desc_index',& if(debug) write(0,*)'in psi_crea_index calling psi_desc_index',&
& size(index_out) & size(index_out)
call psi_desc_index(desc_a%matrix_data,index_in,dep_list(1:,me),& call psi_desc_index(desc_a%matrix_data,index_in,dep_list(1:,me),&
& length_dl(me),desc_a%loc_to_glob,desc_a%glob_to_loc,& & length_dl(me),nsnd,nrcv,desc_a%loc_to_glob,desc_a%glob_to_loc,&
& index_out,glob_idx,info) & index_out,glob_idx,info)
if(debug) write(0,*)'out of psi_desc_index',&
& size(index_out)
nxch = length_dl(me)
if(info /= 0) then if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_desc_index') call psb_errpush(4010,name,a_err='psi_desc_index')
goto 9999 goto 9999

@ -45,58 +45,128 @@ subroutine psi_crea_ovr_elem(desc_overlap,ovr_elem)
integer :: psi_exist_ovr_elem,dim integer :: psi_exist_ovr_elem,dim
external :: psi_exist_ovr_elem external :: psi_exist_ovr_elem
logical, parameter :: usetree=.true. integer :: nel, ip, ix, iel, insize
integer, allocatable :: telem(:,:)
logical, parameter :: usetree=.false.
if (associated(ovr_elem)) then if (associated(ovr_elem)) then
dim_ovr_elem=size(ovr_elem) dim_ovr_elem=size(ovr_elem)
else else
dim_ovr_elem = 0 dim_ovr_elem = 0
endif endif
i=1
pnt_new_elem=1
if (usetree) call initpairsearchtree(pairtree,info) if (usetree) then
do while (desc_overlap(i).ne.-1)
! ...loop over all procs of desc_overlap list.... !
! While running through the column indices exchanged with other procs
i=i+1 ! we have to record them in overlap_elem. We do this by maintaining
do j=1,desc_overlap(i) ! an AVL balanced search tree: at each point counter_e is the next
! ....loop over all overlap indices referred to act proc..... ! free index element. The search routine for gidx will return
if (usetree) then ! glx if gidx was already assigned a local index (glx<counter_e)
! but if gidx was a new index for this process, then it creates
! a new pair (gidx,counter_e), and glx==counter_e. In this case we
! need to record this for the overlap exchange. Otherwise it was
! already there, so we need to record one more parnter in the exchange
!
i=1
pnt_new_elem=1
call initpairsearchtree(pairtree,info)
do while (desc_overlap(i).ne.-1)
! ...loop over all procs of desc_overlap list....
i=i+1
do j=1,desc_overlap(i)
! ....loop over all overlap indices referred to act proc.....
call searchinskeyval(pairtree,desc_overlap(i+j),pnt_new_elem,& call searchinskeyval(pairtree,desc_overlap(i+j),pnt_new_elem,&
& ret,info) & ret,info)
if (ret == pnt_new_elem) ret=-1 if (ret == pnt_new_elem) ret=-1
else if (ret.eq.-1) then
ret=psi_exist_ovr_elem(ovr_elem,pnt_new_elem-2,&
& desc_overlap(i+j)) ! ...this point not exist in ovr_elem list:
endif ! add to it.............................
if (ret.eq.-1) then ! ...check if overflow element_d array......
if ((pnt_new_elem +2) > dim_ovr_elem) then
! ...this point not exist in ovr_elem list: dim_ovr_elem=max(((3*dim_ovr_elem)/2+2),pnt_new_elem+100)
! add to it............................. call psb_realloc(dim_ovr_elem,ovr_elem,info)
ovr_elem(pnt_new_elem)=desc_overlap(i+j) if (info /= 0) then
ovr_elem(pnt_new_elem+1)=2 write(0,*) 'Error in CREA_OVR_ELEM'
pnt_new_elem=pnt_new_elem+2 endif
! ...check if overflow element_d array......
if (pnt_new_elem.gt.dim_ovr_elem) then
dim_ovr_elem=max(((3*size(ovr_elem))/2+2),pnt_new_elem+100)
!!$ write(0,*) 'calling realloc crea_ovr_elem',dim
call psb_realloc(dim_ovr_elem,ovr_elem,info)
if (info /= 0) then
write(0,*) 'Error in CREA_OVR_ELEM'
endif endif
ovr_elem(pnt_new_elem)=desc_overlap(i+j)
ovr_elem(pnt_new_elem+1)=2
pnt_new_elem=pnt_new_elem+2
else
! ....this point already exist in ovr_elem list
! its position is ret............................
ovr_elem(ret+1)=ovr_elem(ret+1)+1
endif endif
else enddo
! ....this point already exist in ovr_elem list i=i+2*desc_overlap(i)+2
! its position is ret............................
ovr_elem(ret+1)=ovr_elem(ret+1)+1
endif
enddo enddo
i=i+2*desc_overlap(i)+2
enddo
! ...add -1 at the end of output list...... ! Add -1 at the end of output list.
ovr_elem(pnt_new_elem)=-1 ! And fix the size to the minimum necessary.
if (usetree) call freepairsearchtree(pairtree) dim_ovr_elem=pnt_new_elem
call psb_realloc(dim_ovr_elem,ovr_elem,info)
if (info /= 0) then
write(0,*) 'Error in CREA_OVR_ELEM'
endif
ovr_elem(pnt_new_elem)=-1
call freepairsearchtree(pairtree)
else
insize = size(desc_overlap)
insize = max(1,(insize+1)/2)
allocate(telem(insize,2),stat=info)
if (info /= 0) then
write(0,*) 'Error in CREA_OVR_ELEM'
endif
i = 1
nel = 0
do while (desc_overlap(i).ne.-1)
! ...loop over all procs of desc_overlap list....
i=i+1
do j=1,desc_overlap(i)
nel = nel + 1
telem(nel,1) = desc_overlap(i+j)
enddo
i=i+2*desc_overlap(i)+2
enddo
if (nel > 0) then
call imsr(nel,telem(:,1))
iel = telem(1,1)
telem(1,2) = 2
ix = 1
ip = 2
do
if (ip > nel) exit
if (telem(ip,1) == iel) then
telem(ix,2) = telem(ix,2) + 1
else
ix = ix + 1
telem(ix,1) = telem(ip,1)
iel = telem(ip,1)
telem(ix,2) = 2
end if
ip = ip + 1
end do
else
ix = 0
end if
dim_ovr_elem=2*ix+1
call psb_realloc(dim_ovr_elem,ovr_elem,info)
iel = 1
do i=1, ix
ovr_elem(iel) = telem(i,1)
ovr_elem(iel+1) = telem(i,2)
iel = iel + 2
end do
ovr_elem(iel) = -1
deallocate(telem)
endif
end subroutine psi_crea_ovr_elem end subroutine psi_crea_ovr_elem

@ -29,7 +29,7 @@
!!$ !!$
!!$ !!$
subroutine psi_desc_index(desc_data,index_in,dep_list,& subroutine psi_desc_index(desc_data,index_in,dep_list,&
& length_dl,loc_to_glob,glob_to_loc,desc_index,& & length_dl,nsnd,nrcv,loc_to_glob,glob_to_loc,desc_index,&
& isglob_in,info) & isglob_in,info)
use psb_realloc_mod use psb_realloc_mod
@ -39,23 +39,22 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
use psb_penv_mod use psb_penv_mod
implicit none implicit none
!c ...array parameters..... ! ...array parameters.....
integer :: desc_data(:),index_in(:),dep_list(:) integer :: desc_data(:),index_in(:),dep_list(:)
integer :: loc_to_glob(:),glob_to_loc(:) integer :: loc_to_glob(:),glob_to_loc(:)
integer,pointer :: desc_index(:) integer,pointer :: desc_index(:)
integer :: length_dl, info integer :: length_dl, nsnd,nrcv,info
logical :: isglob_in logical :: isglob_in
!c ....local scalars... ! ....local scalars...
integer :: j,me,np,i,proc integer :: j,me,np,i,proc
!c ...parameters... ! ...parameters...
integer :: ictxt integer :: ictxt
integer :: no_comm,err integer :: no_comm,err
parameter (no_comm=-1) parameter (no_comm=-1)
!c ...local arrays.. ! ...local arrays..
integer,pointer :: brvindx(:),rvsz(:),& integer,pointer :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:), sndbuf(:), rcvbuf(:)
& bsdindx(:),sdsz(:), sndbuf(:), rcvbuf(:)
integer :: ihinsz,ntot,k,err_act,& integer :: ihinsz,ntot,k,err_act, l_di, &
& idxr, idxs, iszs, iszr, nesd, nerv, icomm & idxr, idxs, iszs, iszr, nesd, nerv, icomm
logical,parameter :: debug=.false., usempi=.true. logical,parameter :: debug=.false., usempi=.true.
@ -65,9 +64,8 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
name='psi_desc_index' name='psi_desc_index'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
!c if mode == 1 then we can use glob_to_loc array
!c else we can't utilize it
ictxt=desc_data(psb_ctxt_) ictxt=desc_data(psb_ctxt_)
icomm=desc_data(psb_mpi_c_)
call psb_info(ictxt,me,np) call psb_info(ictxt,me,np)
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
@ -80,13 +78,12 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif endif
call psb_get_mpicomm(ictxt,icomm) !
!c ! first, find out the total sizes to be exchanged.
!c first, find out the total sizes to be exchanged. ! note: things marked here as sndbuf/rcvbuf (for mpi) corresponds to things
!c note: things marked here as sndbuf/rcvbuf (for mpi) corresponds to things ! to be received/sent (in the final psblas descriptor).
!c to be received/sent (in the final psblas descriptor). ! be careful of the inversion
!c be careful of the inversion !
!c
allocate(sdsz(np),rvsz(np),bsdindx(np),brvindx(np),stat=info) allocate(sdsz(np),rvsz(np),bsdindx(np),brvindx(np),stat=info)
if(info /= 0) then if(info /= 0) then
info=4000 info=4000
@ -126,6 +123,9 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
end do end do
iszs = sum(sdsz) iszs = sum(sdsz)
iszr = sum(rvsz) iszr = sum(rvsz)
nsnd = iszr
nrcv = iszs
if ((iszs /= idxs).or.(iszr /= idxr)) then if ((iszs /= idxs).or.(iszr /= idxr)) then
write(0,*) 'strange results???', iszs,idxs,iszr,idxr write(0,*) 'strange results???', iszs,idxs,iszr,idxr
end if end if
@ -134,11 +134,14 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif endif
ntot = (3*(max(count(sdsz>0),count(rvsz>0)))+ iszs + iszr) + 1 ntot = (3*(count((sdsz>0).or.(rvsz>0)))+ iszs + iszr) + 1
if (size(desc_index) < ntot) then
!c$$$ write(0,*) 'potential error on desc_index :', if (associated(desc_index)) then
!c$$$ + length_dh, size(desc_index),ntot l_di = size(desc_index)
write(0,*) 'calling irealloc psi_desc_index ',ntot else
l_di = 0
endif
if (l_di < ntot) then
call psb_realloc(ntot,desc_index,info) call psb_realloc(ntot,desc_index,info)
endif endif
if (info /= 0) then if (info /= 0) then
@ -167,9 +170,9 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
proc = index_in(i) proc = index_in(i)
i = i + 1 i = i + 1
nerv = index_in(i) nerv = index_in(i)
! c !
! c note that here bsdinx is zero-based, hence the following loop ! note that here bsdinx is zero-based, hence the following loop
! c !
if (isglob_in) then if (isglob_in) then
do j=1, nerv do j=1, nerv
sndbuf(bsdindx(proc+1)+j) = (index_in(i+j)) sndbuf(bsdindx(proc+1)+j) = (index_in(i+j))
@ -187,9 +190,9 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
write(0,*) me,' prepared send buffer ' write(0,*) me,' prepared send buffer '
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif endif
!c !
!c now have to regenerate bsdindx ! now have to regenerate bsdindx
!c !
idxs = 0 idxs = 0
idxr = 0 idxr = 0
do i=1, length_dl do i=1, length_dl
@ -207,10 +210,10 @@ subroutine psi_desc_index(desc_data,index_in,dep_list,&
goto 9999 goto 9999
end if end if
!c !
!c at this point we can finally build the output desc_index. beware ! at this point we can finally build the output desc_index. beware
!c of snd/rcv inversion. ! of snd/rcv inversion.
!c !
i = 1 i = 1
do k = 1, length_dl do k = 1, length_dl
proc = dep_list(k) proc = dep_list(k)

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

@ -119,8 +119,8 @@ c desc_str list.
c length_dl integer array(0:np) c length_dl integer array(0:np)
c length_dl(i) is the length of dep_list(*,i) list c length_dl(i) is the length of dep_list(*,i) list
use psb_penv_mod use psb_penv_mod
use psb_const_mod
implicit none implicit none
include 'psb_const.fh'
include 'mpif.h' include 'mpif.h'
c ....scalar parameters... c ....scalar parameters...
integer np,dl_lda,mode, info integer np,dl_lda,mode, info

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

@ -19,20 +19,22 @@ LIBDIR = ../../lib
psb_realloc_mod.o : psb_error_mod.o psb_realloc_mod.o : psb_error_mod.o
psb_spmat_type.o : psb_realloc_mod.o psb_const_mod.o psb_spmat_type.o : psb_realloc_mod.o psb_const_mod.o
psb_error_mod.o: psb_const_mod.o psb_error_mod.o: psb_const_mod.o
psb_const_mod.f90: psb_const.fh psb_penv_mod.o: psb_const_mod.o psb_error_mod.o
psb_penv_mod.o : psb_const_mod.o psb_error_mod.o psi_mod.o: psb_penv_mod.o psb_error_mod.o psb_desc_type.o
psb_methd_mod.o: psb_serial_mod.o psb_desc_type.o psb_prec_type.o
psb_sparse_mod.o: $(MODULES) $(MPFOBJS)
lib: $(MODULES) mpfobjs $(OBJS) lib: mpfobjs $(MODULES) $(OBJS)
$(AR) $(LIBDIR)/$(LIBNAME) $(MODULES) $(OBJS) $(MPFOBJS) $(AR) $(LIBDIR)/$(LIBNAME) $(MODULES) $(OBJS) $(MPFOBJS)
$(RANLIB) $(LIBDIR)/$(LIBNAME) $(RANLIB) $(LIBDIR)/$(LIBNAME)
cp *$(.mod) ./psb_const.fh ./parts.fh ../../lib /bin/cp *$(.mod) ./parts.fh ../../lib
mpfobjs: mpfobjs:
(make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)") (make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)")
clean: clean:
/bin/rm -f $(MODULES) $(OBJS) *$(.mod) /bin/rm -f $(MODULES) $(OBJS) $(MPFOBJS) *$(.mod)
veryclean: clean veryclean: clean

@ -57,16 +57,24 @@
integer, parameter :: psb_dec_type_=1, psb_m_=2,psb_n_=3 integer, parameter :: psb_dec_type_=1, psb_m_=2,psb_n_=3
integer, parameter :: psb_n_row_=4, psb_n_col_=5,psb_ctxt_=6 integer, parameter :: psb_n_row_=4, psb_n_col_=5,psb_ctxt_=6
integer, parameter :: psb_loc_to_glob_=7 integer, parameter :: psb_loc_to_glob_=7
integer, parameter :: psb_mpi_c_=9,psb_mdata_size_=10 integer, parameter :: psb_thal_xch_=11
integer, parameter :: psb_thal_snd_=12
integer, parameter :: psb_thal_rcv_=13
integer, parameter :: psb_tovr_xch_=14
integer, parameter :: psb_tovr_snd_=15
integer, parameter :: psb_tovr_rcv_=16
integer, parameter :: psb_mpi_c_=9,psb_mdata_size_=20
integer, parameter :: psb_desc_asb_=3099 integer, parameter :: psb_desc_asb_=3099
integer, parameter :: psb_desc_bld_=psb_desc_asb_+1 integer, parameter :: psb_desc_bld_=psb_desc_asb_+1
integer, parameter :: psb_desc_repl_=3199 integer, parameter :: psb_desc_repl_=3199
integer, parameter :: psb_desc_upd_=psb_desc_bld_+1 integer, parameter :: psb_desc_upd_=psb_desc_bld_+1
integer, parameter :: psb_desc_upd_asb_=psb_desc_upd_+1 integer, parameter :: psb_desc_upd_asb_=psb_desc_upd_+1
integer, parameter :: psb_hpnt_nentries_=7
! !
! Constants for desc_a handling ! Constants for desc_a handling
! !
integer, parameter :: psb_upd_glbnum_=998 integer, parameter :: psb_upd_glbnum_=998
integer, parameter :: psb_upd_locnum_=997 integer, parameter :: psb_upd_locnum_=997
integer, parameter :: psb_proc_id_=0, psb_n_elem_recv_=1 integer, parameter :: psb_proc_id_=0, psb_n_elem_recv_=1

@ -31,6 +31,109 @@
module psb_const_mod module psb_const_mod
include 'psb_const.fh' !
! Communication, prolongation & restriction
end module psb_const_mod !
integer, parameter :: psb_nohalo_=0, psb_halo_=4
integer, parameter :: psb_none_=0, psb_sum_=1
integer, parameter :: psb_avg_=2, psb_square_root_=3
integer, parameter :: psb_swap_send_=1, psb_swap_recv_=2
integer, parameter :: psb_swap_sync_=4, psb_swap_mpi_=8
!
! Data checks
!
integer, parameter :: psb_deadlock_check_=0
integer, parameter :: psb_local_mtrx_check_=1
integer, parameter :: psb_local_comm_check_=2
integer, parameter :: psb_consistency_check_=3
integer, parameter :: psb_global_check_=4
integer, parameter :: psb_order_communication_=5
integer, parameter :: psb_change_represent_=6
integer, parameter :: psb_loc_to_glob_check_=7
integer, parameter :: psb_convert_halo_=1, psb_convert_ovrlap_=2
integer, parameter :: psb_act_ret_=0, psb_act_abort_=1, no_err_=0
!
! Entries and values in desc%matrix_data
!
integer, parameter :: psb_dec_type_=1, psb_m_=2,psb_n_=3
integer, parameter :: psb_n_row_=4, psb_n_col_=5,psb_ctxt_=6
integer, parameter :: psb_loc_to_glob_=7
integer, parameter :: psb_thal_xch_=11
integer, parameter :: psb_thal_snd_=12
integer, parameter :: psb_thal_rcv_=13
integer, parameter :: psb_tovr_xch_=14
integer, parameter :: psb_tovr_snd_=15
integer, parameter :: psb_tovr_rcv_=16
integer, parameter :: psb_mpi_c_=9,psb_mdata_size_=20
integer, parameter :: psb_desc_asb_=3099
integer, parameter :: psb_desc_bld_=psb_desc_asb_+1
integer, parameter :: psb_desc_repl_=3199
integer, parameter :: psb_desc_upd_=psb_desc_bld_+1
integer, parameter :: psb_desc_upd_asb_=psb_desc_upd_+1
integer, parameter :: psb_hpnt_nentries_=7
!
! Constants for desc_a handling
!
integer, parameter :: psb_upd_glbnum_=998
integer, parameter :: psb_upd_locnum_=997
integer, parameter :: psb_proc_id_=0, psb_n_elem_recv_=1
integer, parameter :: psb_elem_recv_=2, psb_n_elem_send_=2
integer, parameter :: psb_elem_send_=3, psb_n_ovrlp_elem_=1
integer, parameter :: psb_ovrlp_elem_to_=2, psb_ovrlp_elem_=0
integer, parameter :: psb_n_dom_ovr_=1
integer, parameter :: psb_no_comm_=-1
integer, parameter :: psb_comm_halo_=0, psb_comm_ovr_=1
!
! Queries into spmat%info
!
integer, parameter :: psb_nztotreq_=1, psb_nzrowreq_=2
integer, parameter :: psb_nzsizereq_=3
!
! Entries and values for spmat%info
!
integer, parameter :: psb_nnz_=1
integer, parameter :: psb_del_bnd_=7, psb_srtd_=8
integer, parameter :: psb_state_=9
integer, parameter :: psb_upd_pnt_=10
integer, parameter :: psb_dupl_=11, psb_upd_=12
integer, parameter :: psb_ifasize_=16
integer, parameter :: psb_spmat_null_=0, psb_spmat_bld_=1
integer, parameter :: psb_spmat_asb_=2, psb_spmat_upd_=4
integer, parameter :: psb_ireg_flgs_=10, psb_ip2_=0
integer, parameter :: psb_iflag_=2, psb_ichk_=3
integer, parameter :: psb_nnzt_=4, psb_zero_=5,psb_ipc_=6
integer, parameter :: psb_dupl_ovwrt_ = 0
integer, parameter :: psb_dupl_add_ = 1
integer, parameter :: psb_dupl_err_ = 2
integer, parameter :: psb_dupl_def_ = psb_dupl_ovwrt_
integer, parameter :: psb_upd_dflt_ = 0
integer, parameter :: psb_upd_srch_ = 98764
integer, parameter :: psb_upd_perm_ = 98765
integer, parameter :: psb_isrtdcoo_ = 98761
integer, parameter :: psb_maxjdrows_=8, psb_minjdrows_=4
integer, parameter :: psb_dbleint_=2
!
! Error handling
!
integer, parameter :: act_ret=0, act_abort=1, no_err=0
!
! Handy & miscellaneous constants
!
integer, parameter :: ione=1, izero=0
integer, parameter :: itwo=2, ithree=3,mone=-1, psb_root_=0
real(kind(1.d0)), parameter :: psb_colrow_=0.33, psb_percent_=0.7
real(kind(1.d0)), parameter :: dzero=0.d0, done=1.d0
complex(kind(1.d0)), parameter :: zzero=(0.d0,0.0d0)
complex(kind(1.d0)), parameter :: zone=(1.d0,0.0d0)
real(kind(1.d0)), parameter :: epstol=1.d-32
character, parameter :: psb_all_='A', psb_topdef_=' '
character(len=5) :: psb_fidef_='CSR'
end module psb_const_mod

@ -42,45 +42,47 @@ module psb_descriptor_type
! desc_type contains data for communications. ! desc_type contains data for communications.
type psb_desc_type type psb_desc_type
! contain decomposition informations ! contain decomposition informations
integer, pointer :: matrix_data(:)=>null() integer, pointer :: matrix_data(:)=>null()
! contain index of halo elements to send/receive ! contain index of halo elements to send/receive
integer, pointer :: halo_index(:)=>null() integer, pointer :: halo_index(:)=>null()
! contain indices of boundary elements ! contain indices of boundary elements
integer, pointer :: bnd_elem(:)=>null() integer, pointer :: bnd_elem(:)=>null()
! contain index of overlap elements to send/receive ! contain index of overlap elements to send/receive
integer, pointer :: ovrlap_elem(:)=>null() integer, pointer :: ovrlap_elem(:)=>null()
! contain for each local overlap element, the number of times ! contain for each local overlap element, the number of times
! that is duplicated ! that is duplicated
integer, pointer :: ovrlap_index(:)=>null() integer, pointer :: ovrlap_index(:)=>null()
! contain for each local element the corresponding global index ! contain for each local element the corresponding global index
integer, pointer :: loc_to_glob(:)=>null() integer, pointer :: loc_to_glob(:)=>null()
! contain for each global element the corresponding local index, ! contain for each global element the corresponding local index,
! if exist. ! if exist.
integer, pointer :: glob_to_loc (:)=>null() integer, pointer :: glob_to_loc (:)=>null()
! local renumbering induced by sparse matrix storage. ! local renumbering induced by sparse matrix storage.
integer, pointer :: lprm(:)=>null() integer, pointer :: lprm(:)=>null()
! index space in case it is not just the contiguous range 1:n ! index space in case it is not just the contiguous range 1:n
integer, pointer :: idx_space(:)=>null() integer, pointer :: idx_space(:)=>null()
end type psb_desc_type end type psb_desc_type
contains contains
subroutine psb_nullify_desc(desc) subroutine psb_nullify_desc(desc)
type(psb_desc_type), intent(inout) :: desc type(psb_desc_type), intent(inout) :: desc
nullify(desc%matrix_data,desc%loc_to_glob,desc%glob_to_loc,& nullify(desc%matrix_data,desc%loc_to_glob,desc%glob_to_loc,&
&desc%halo_index,desc%bnd_elem,desc%ovrlap_elem,& & desc%halo_index,desc%bnd_elem,desc%ovrlap_elem,&
&desc%ovrlap_index, desc%lprm, desc%idx_space) & desc%ovrlap_index, desc%lprm, desc%idx_space)!,&
! & desc%halo_pt,desc%ovrlap_pt)
end subroutine psb_nullify_desc end subroutine psb_nullify_desc
logical function psb_is_ok_dec(dectype) logical function psb_is_ok_dec(dectype)
integer :: dectype integer :: dectype
psb_is_ok_dec = ((dectype == psb_desc_asb_).or.(dectype == psb_desc_bld_).or.& psb_is_ok_dec = ((dectype == psb_desc_asb_).or.(dectype == psb_desc_bld_).or.&
&(dectype == psb_desc_upd_).or.(dectype== psb_desc_upd_asb_)) & (dectype == psb_desc_upd_).or.(dectype== psb_desc_upd_asb_).or.&
& (dectype== psb_desc_repl_))
end function psb_is_ok_dec end function psb_is_ok_dec
logical function psb_is_bld_dec(dectype) logical function psb_is_bld_dec(dectype)
@ -106,7 +108,8 @@ contains
logical function psb_is_asb_dec(dectype) logical function psb_is_asb_dec(dectype)
integer :: dectype integer :: dectype
psb_is_asb_dec = (dectype == psb_desc_asb_) psb_is_asb_dec = (dectype == psb_desc_asb_).or.&
& (dectype== psb_desc_repl_)
end function psb_is_asb_dec end function psb_is_asb_dec

@ -352,9 +352,9 @@ contains
case(150) case(150)
write (0,'("indices in input array are not belonging to the calling process ",i0)')i_e_d(1) write (0,'("indices in input array are not belonging to the calling process ",i0)')i_e_d(1)
case(290) case(290)
write (0,'("Is not possible to call this routine without calling before psdalloc on the same matrix")') write (0,'("To call this routine you must first call psb_geall on the same matrix")')
case(295) case(295)
write (0,'("Is not possible to call this routine without calling before psdspalloc on the same matrix")') write (0,'("To call this routine you must first call psb_spall on the same matrix")')
case(300) case(300)
write (0,'("Input argument n. ",i0," must be equal to entry n. ",i0," in array input argument n.",i0)') & write (0,'("Input argument n. ",i0," must be equal to entry n. ",i0," in array input argument n.",i0)') &
& i_e_d(1),i_e_d(4),i_e_d(3) & i_e_d(1),i_e_d(4),i_e_d(3)
@ -376,7 +376,7 @@ contains
case(600) case(600)
write (0,'("Sparse Matrix and decsriptors are in an invalid state for this subroutine call: ",i0)')i_e_d(1) write (0,'("Sparse Matrix and decsriptors are in an invalid state for this subroutine call: ",i0)')i_e_d(1)
case (1122) case (1122)
write (0,'("Invalid state for DESC_A")') write (0,'("Invalid state for communication descriptor")')
case (1123) case (1123)
write (0,'("Invalid combined state for A and DESC_A")') write (0,'("Invalid combined state for A and DESC_A")')
case(1124:1999) case(1124:1999)

@ -44,7 +44,287 @@ module psb_realloc_mod
module procedure psb_dreallocatez2 module procedure psb_dreallocatez2
end Interface end Interface
Contains Interface psb_safe_cpy
module procedure psb_icpy1d,psb_icpy2d, &
& psb_dcpy1d, psb_dcpy2d, psb_zcpy1d, psb_zcpy2d
end Interface
contains
subroutine psb_icpy1d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
Integer,pointer :: vin(:),vout(:)
integer :: info
! ...Local Variables
Integer :: isz,err_act
character(len=20) :: name, char_err
logical, parameter :: debug=.false.
name='psb_cpy1d'
call psb_erractionsave(err_act)
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
isz = size(vin)
call psb_realloc(isz,vout,info)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
vout(:) = vin(:)
endif
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
else
call psb_error()
end if
return
end subroutine psb_icpy1d
subroutine psb_icpy2d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
Integer,pointer :: vin(:,:),vout(:,:)
integer :: info
! ...Local Variables
Integer :: isz1, isz2,err_act
character(len=20) :: name, char_err
logical, parameter :: debug=.false.
name='psb_cpy1d'
call psb_erractionsave(err_act)
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
isz1 = size(vin,1)
isz2 = size(vin,2)
call psb_realloc(isz1,isz2,vout,info)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
vout(:,:) = vin(:,:)
endif
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
else
call psb_error()
end if
return
end subroutine psb_icpy2d
subroutine psb_dcpy1d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
real(kind(1.d0)),pointer :: vin(:),vout(:)
integer :: info
! ...Local Variables
Integer :: isz,err_act
character(len=20) :: name, char_err
logical, parameter :: debug=.false.
name='psb_cpy1d'
call psb_erractionsave(err_act)
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
isz = size(vin)
call psb_realloc(isz,vout,info)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
vout(:) = vin(:)
endif
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
else
call psb_error()
end if
return
end subroutine psb_dcpy1d
subroutine psb_dcpy2d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
real(kind(1.d0)),pointer :: vin(:,:),vout(:,:)
integer :: info
! ...Local Variables
Integer :: isz1, isz2,err_act
character(len=20) :: name, char_err
logical, parameter :: debug=.false.
name='psb_cpy1d'
call psb_erractionsave(err_act)
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
isz1 = size(vin,1)
isz2 = size(vin,2)
call psb_realloc(isz1,isz2,vout,info)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
vout(:,:) = vin(:,:)
endif
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
else
call psb_error()
end if
return
end subroutine psb_dcpy2d
subroutine psb_zcpy1d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
complex(kind(1.d0)),pointer :: vin(:),vout(:)
integer :: info
! ...Local Variables
Integer :: isz,err_act
character(len=20) :: name, char_err
logical, parameter :: debug=.false.
name='psb_cpy1d'
call psb_erractionsave(err_act)
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
isz = size(vin)
call psb_realloc(isz,vout,info)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
vout(:) = vin(:)
endif
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
else
call psb_error()
end if
return
end subroutine psb_zcpy1d
subroutine psb_zcpy2d(vin,vout,info)
use psb_error_mod
! ...Subroutine Arguments
complex(kind(1.d0)),pointer :: vin(:,:),vout(:,:)
integer :: info
! ...Local Variables
Integer :: isz1, isz2,err_act
character(len=20) :: name, char_err
logical, parameter :: debug=.false.
name='psb_cpy1d'
call psb_erractionsave(err_act)
if(psb_get_errstatus().ne.0) return
info = 0
if (associated(vin)) then
isz1 = size(vin,1)
isz2 = size(vin,2)
call psb_realloc(isz1,isz2,vout,info)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
vout(:,:) = vin(:,:)
endif
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
else
call psb_error()
end if
return
end subroutine psb_zcpy2d
Subroutine psb_dreallocate1i(len,rrax,info,pad) Subroutine psb_dreallocate1i(len,rrax,info,pad)
use psb_error_mod use psb_error_mod
@ -120,9 +400,6 @@ Contains
End Subroutine psb_dreallocate1i End Subroutine psb_dreallocate1i
Subroutine psb_dreallocate1d(len,rrax,info,pad) Subroutine psb_dreallocate1d(len,rrax,info,pad)
use psb_error_mod use psb_error_mod
@ -154,7 +431,6 @@ Contains
goto 9999 goto 9999
end if end if
m = min(dim,len) m = min(dim,len)
!!$ write(0,*) 'DA: copying ',min(len,dim)
if (.true.) then if (.true.) then
do i=1,m do i=1,m
tmp(i) = rrax(i) tmp(i) = rrax(i)
@ -162,7 +438,6 @@ Contains
else else
tmp(1:m) = rrax(1:m) tmp(1:m) = rrax(1:m)
end if end if
!!$ write(0,*) 'DA: copying done ',m
Deallocate(rrax,stat=info) Deallocate(rrax,stat=info)
if (info /= 0) then if (info /= 0) then
err=4000 err=4000

@ -393,27 +393,18 @@ contains
Integer, intent(out) :: info Integer, intent(out) :: info
!locals !locals
Integer :: nza,nz1, nz2, nzl, nzr
logical, parameter :: debug=.false.
INFO = 0 INFO = 0
call psb_nullify_sp(b)
nza = size(a%aspk) call psb_safe_cpy(a%aspk,b%aspk,info)
nz1 = size(a%ia1) if (info == 0) call psb_safe_cpy(a%ia1,b%ia1,info)
nz2 = size(a%ia2) if (info == 0) call psb_safe_cpy(a%ia2,b%ia2,info)
nzl = size(a%pl) if (info == 0) call psb_safe_cpy(a%pl,b%pl,info)
nzr = size(a%pr) if (info == 0) call psb_safe_cpy(a%pr,b%pr,info)
allocate(b%aspk(nza),b%ia1(nz1),b%ia2(nz2),&
& b%pl(nzl),b%pr(nzr),stat=info)
if (info /= 0) then if (info /= 0) then
info=2023 info=2023
return return
Endif Endif
b%aspk(:) = a%aspk(:)
b%ia1(:) = a%ia1(:)
b%ia2(:) = a%ia2(:)
b%pl(:) = a%pl(:)
b%pr(:) = a%pr(:)
b%infoa(:) = a%infoa(:) b%infoa(:) = a%infoa(:)
b%fida = a%fida b%fida = a%fida
b%descra = a%descra b%descra = a%descra
@ -913,27 +904,21 @@ contains
Integer, intent(out) :: info Integer, intent(out) :: info
!locals !locals
Integer :: nza,nz1, nz2, nzl, nzr
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
INFO = 0
nza = size(a%aspk) INFO = 0
nz1 = size(a%ia1) call psb_nullify_sp(b)
nz2 = size(a%ia2) call psb_safe_cpy(a%aspk,b%aspk,info)
nzl = size(a%pl) if (info == 0) call psb_safe_cpy(a%ia1,b%ia1,info)
nzr = size(a%pr) if (info == 0) call psb_safe_cpy(a%ia2,b%ia2,info)
allocate(b%aspk(nza),b%ia1(nz1),b%ia2(nz2),& if (info == 0) call psb_safe_cpy(a%pl,b%pl,info)
& b%pl(nzl),b%pr(nzr),stat=info) if (info == 0) call psb_safe_cpy(a%pr,b%pr,info)
if (info /= 0) then if (info /= 0) then
info=2023 info=2023
return return
Endif Endif
b%aspk(:) = a%aspk(:)
b%ia1(:) = a%ia1(:)
b%ia2(:) = a%ia2(:)
b%pl(:) = a%pl(:)
b%pr(:) = a%pr(:)
b%infoa(:) = a%infoa(:) b%infoa(:) = a%infoa(:)
b%fida = a%fida b%fida = a%fida
b%descra = a%descra b%descra = a%descra

@ -51,10 +51,10 @@ module psi_mod
end interface end interface
interface interface
subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,info) subroutine psi_crea_index(desc_a,index_in,index_out,glob_idx,nxch,nsnd,nrcv,info)
use psb_descriptor_type use psb_descriptor_type
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info integer, intent(out) :: info,nxch,nsnd,nrcv
integer, intent(in) :: index_in(:) integer, intent(in) :: index_in(:)
integer, pointer :: index_out(:) integer, pointer :: index_out(:)
logical :: glob_idx logical :: glob_idx
@ -70,12 +70,12 @@ module psi_mod
interface interface
subroutine psi_desc_index(desc_data,index_in,dep_list,& subroutine psi_desc_index(desc_data,index_in,dep_list,&
& length_dl,loc_to_glob,glob_to_loc,desc_index,& & length_dl,nsnd,nrcv,loc_to_glob,glob_to_loc,desc_index,&
& isglob_in,info) & isglob_in,info)
integer :: desc_data(:),index_in(:),dep_list(:) integer :: desc_data(:),index_in(:),dep_list(:)
integer :: loc_to_glob(:),glob_to_loc(:) integer :: loc_to_glob(:),glob_to_loc(:)
integer,pointer :: desc_index(:) integer,pointer :: desc_index(:)
integer :: length_dl, info integer :: length_dl,nsnd,nrcv, info
logical :: isglob_in logical :: isglob_in
end subroutine psi_desc_index end subroutine psi_desc_index
end interface end interface
@ -238,4 +238,104 @@ module psi_mod
end subroutine psi_zsctv end subroutine psi_zsctv
end interface end interface
interface psi_cnv_dsc
module procedure psi_cnv_dsc
end interface
contains
subroutine psi_cnv_dsc(halo_in,ovrlap_in,cdesc, info)
use psb_const_mod
use psb_error_mod
use psb_penv_mod
use psb_descriptor_type
implicit none
! ....scalars parameters....
integer, intent(in) :: halo_in(:), ovrlap_in(:)
type(psb_desc_type), intent(inout) :: cdesc
integer, intent(out) :: info
! ....local scalars....
integer :: i,np,me,proc, max_index
integer :: ictxt, err_act,nxch,nsnd,nrcv
! ...local array...
integer :: int_err(5)
integer, pointer :: idx_out(:)
! ...parameters
logical, parameter :: debug=.false.
character(len=20) :: name
name='psi_bld_cdesc'
call psb_get_erraction(err_act)
info = 0
ictxt = cdesc%matrix_data(psb_ctxt_)
call psb_info(ictxt,me,np)
if (np == -1) then
info = 2010
call psb_errpush(info,name)
goto 9999
endif
! first the halo index
if (debug) write(0,*) me,'Calling crea_index on halo'
idx_out => null()
call psi_crea_index(cdesc,halo_in, idx_out,.false.,nxch,nsnd,nrcv,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
cdesc%halo_index => idx_out
cdesc%matrix_data(psb_thal_xch_) = nxch
cdesc%matrix_data(psb_thal_snd_) = nsnd
cdesc%matrix_data(psb_thal_rcv_) = nrcv
if (debug) write(0,*) me,'Done crea_index on halo'
if (debug) write(0,*) me,'Calling crea_index on ovrlap'
! then the overlap index
idx_out => null()
call psi_crea_index(cdesc,ovrlap_in, idx_out,.true.,nxch,nsnd,nrcv,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
cdesc%ovrlap_index => idx_out
cdesc%matrix_data(psb_tovr_xch_) = nxch
cdesc%matrix_data(psb_tovr_snd_) = nsnd
cdesc%matrix_data(psb_tovr_rcv_) = nrcv
if (debug) write(0,*) me,'Calling crea_ovr_elem'
! next ovrlap_elem
call psi_crea_ovr_elem(cdesc%ovrlap_index,cdesc%ovrlap_elem)
if (debug) write(0,*) me,'Done crea_ovr_elem'
! finally bnd_elem
call psi_crea_bnd_elem(cdesc,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_bnd_elem')
goto 9999
end if
if (debug) write(0,*) me,'Done crea_bnd_elem'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == act_abort) then
call psb_error(ictxt)
return
end if
return
end subroutine psi_cnv_dsc
end module psi_mod end module psi_mod

@ -243,6 +243,10 @@ contains
if (p%iprcparm(coarse_mat_) == mat_repl_) then if (p%iprcparm(coarse_mat_) == mat_repl_) then
call psb_cdrep(ntaggr,ictxt,desc_p,info) call psb_cdrep(ntaggr,ictxt,desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
nzbr(:) = 0 nzbr(:) = 0
nzbr(me+1) = irs nzbr(me+1) = irs
@ -287,6 +291,12 @@ contains
else if (p%iprcparm(coarse_mat_) == mat_distr_) then else if (p%iprcparm(coarse_mat_) == mat_distr_) then
call psb_cddec(naggr,ictxt,desc_p,info) call psb_cddec(naggr,ictxt,desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cddec')
goto 9999
end if
call psb_sp_clone(b,bg,info) call psb_sp_clone(b,bg,info)
if(info /= 0) then if(info /= 0) then
call psb_errpush(4010,name,a_err='spclone') call psb_errpush(4010,name,a_err='spclone')
@ -916,6 +926,10 @@ contains
nzbr(me+1) = b%infoa(psb_nnz_) nzbr(me+1) = b%infoa(psb_nnz_)
call psb_cdrep(ntaggr,ictxt,desc_p,info) call psb_cdrep(ntaggr,ictxt,desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
call psb_sum(ictxt,nzbr(1:np)) call psb_sum(ictxt,nzbr(1:np))
nzbg = sum(nzbr) nzbg = sum(nzbr)
@ -970,6 +984,10 @@ contains
goto 9999 goto 9999
end if end if
call psb_cddec(naggr,ictxt,desc_p,info) call psb_cddec(naggr,ictxt,desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cddec')
goto 9999
end if
call psb_sp_free(b,info) call psb_sp_free(b,info)
if(info /= 0) then if(info /= 0) then
@ -985,6 +1003,10 @@ contains
nzbr(me+1) = b%infoa(psb_nnz_) nzbr(me+1) = b%infoa(psb_nnz_)
call psb_cdrep(ntaggr,ictxt,desc_p,info) call psb_cdrep(ntaggr,ictxt,desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdrep')
goto 9999
end if
call psb_sum(ictxt,nzbr(1:np)) call psb_sum(ictxt,nzbr(1:np))

@ -31,8 +31,8 @@ C
SUBROUTINE CHECK_DIM(M, N, IA, NG, IA2, SUBROUTINE CHECK_DIM(M, N, IA, NG, IA2,
+ NZ, LARN, LIAN1, LIAN2, IERRV) + NZ, LARN, LIAN1, LIAN2, IERRV)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -34,8 +34,8 @@ c
* p1,descrn,arn,ia1n,ia2n,infon,p2,larn,lia1n, * p1,descrn,arn,ia1n,ia2n,infon,p2,larn,lia1n,
* lia2n,aux,laux,ierror) * lia2n,aux,laux,ierror)
use psb_const_mod
implicit none implicit none
include 'psb_const.fh'
c .. scalar arguments .. c .. scalar arguments ..
integer larn, laux, lia1n, lia2n, integer larn, laux, lia1n, lia2n,

@ -35,8 +35,8 @@ C
* P1,DESCRN,ARN,IAN1,IAN2,INFON,P2,LARN,LIAN1, * P1,DESCRN,ARN,IAN1,IAN2,INFON,P2,LARN,LIAN1,
* LIAN2,AUX,LAUX,IERROR) * LIAN2,AUX,LAUX,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -32,8 +32,8 @@ C
* IP1,DESCRN,ARN,IAN1,IAN2,INFON,IP2,LARN,LIAN1, * IP1,DESCRN,ARN,IAN1,IAN2,INFON,IP2,LARN,LIAN1,
* LIAN2,AUX,LAUX,IERROR) * LIAN2,AUX,LAUX,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -59,8 +59,8 @@ C DESCRN Output matrix in JAD format
C ARN,IAN1 C ARN,IAN1
C IAN2,INFON, IP1, IP2 C IAN2,INFON, IP1, IP2
C C
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -116,8 +116,8 @@ C
C C
SUBROUTINE DCSRP1(TRANS,M,N,DESCRA,JA,IA, SUBROUTINE DCSRP1(TRANS,M,N,DESCRA,JA,IA,
+ P,WORK,IWORK,LWORK,IERROR) + P,WORK,IWORK,LWORK,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C .. Scalar Arguments .. C .. Scalar Arguments ..
INTEGER LWORK,M, N, IERROR INTEGER LWORK,M, N, IERROR
CHARACTER TRANS CHARACTER TRANS

@ -31,8 +31,8 @@ C
SUBROUTINE DGIND_TRI(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA, SUBROUTINE DGIND_TRI(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA,
+ LARN,LKA,LJA,IPERM,WORK, LWORK, IERROR) + LARN,LKA,LJA,IPERM,WORK, LWORK, IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C ... Scalar arguments ... C ... Scalar arguments ...

@ -31,8 +31,8 @@ C
SUBROUTINE DGINDEX(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA, INFON, SUBROUTINE DGINDEX(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA, INFON,
+ LARN,LKA,LJA,IPERM,WORK, LWORK, SIZE_REQ, IERROR) + LARN,LKA,LJA,IPERM,WORK, LWORK, SIZE_REQ, IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C ... Scalar arguments ... C ... Scalar arguments ...
INTEGER M, LWORK,N_BLOCKS,LARN,LKA,LJA, INTEGER M, LWORK,N_BLOCKS,LARN,LKA,LJA,

@ -31,8 +31,8 @@ C
SUBROUTINE DJDCO(TRANS,M,N,DESCRA,AR,IA1,IA2,IPERM,INFO, SUBROUTINE DJDCO(TRANS,M,N,DESCRA,AR,IA1,IA2,IPERM,INFO,
* IP1,DESCRN,ARN,IA1N,IA2N,INFON,IP2,LARN,LIA1N, * IP1,DESCRN,ARN,IA1N,IA2N,INFON,IP2,LARN,LIA1N,
* LIA2N,AUX,LAUX,IERROR) * LIA2N,AUX,LAUX,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..
INTEGER LARN, LAUX, LIA1N, LIA2N, M, N, IERROR INTEGER LARN, LAUX, LIA1N, LIA2N, M, N, IERROR

@ -36,8 +36,8 @@ C
* IP1,DESCRN,ARN,IA1N,IA2N,INFON,IP2,LARN,LIA1N, * IP1,DESCRN,ARN,IA1N,IA2N,INFON,IP2,LARN,LIA1N,
* LIA2N,AUX,LAUX,IERROR) * LIA2N,AUX,LAUX,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -29,9 +29,9 @@ C POSSIBILITY OF SUCH DAMAGE.
C C
C C
SUBROUTINE GEN_BLOCK(M,NG,IA,AUX) SUBROUTINE GEN_BLOCK(M,NG,IA,AUX)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
INTEGER M, NG INTEGER M, NG
INTEGER IA(3,*), AUX(*) INTEGER IA(3,*), AUX(*)

@ -29,9 +29,9 @@ C POSSIBILITY OF SUCH DAMAGE.
C C
C C
SUBROUTINE PARTITION(M, WORK, IA, N_BLOCK) SUBROUTINE PARTITION(M, WORK, IA, N_BLOCK)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C ...Scalar arguments... C ...Scalar arguments...

@ -34,8 +34,8 @@ c
* p1,descrn,arn,ia1n,ia2n,infon,p2,larn,lia1n, * p1,descrn,arn,ia1n,ia2n,infon,p2,larn,lia1n,
* lia2n,aux,laux,ierror) * lia2n,aux,laux,ierror)
use psb_const_mod
implicit none implicit none
include 'psb_const.fh'
c .. scalar arguments .. c .. scalar arguments ..
integer larn, laux, lia1n, lia2n, integer larn, laux, lia1n, lia2n,

@ -35,8 +35,8 @@ C
* P1,DESCRN,ARN,IAN1,IAN2,INFON,P2,LARN,LIAN1, * P1,DESCRN,ARN,IAN1,IAN2,INFON,P2,LARN,LIAN1,
* LIAN2,AUX,LAUX,IERROR) * LIAN2,AUX,LAUX,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -32,8 +32,8 @@ C
* IP1,DESCRN,ARN,IAN1,IAN2,INFON,IP2,LARN,LIAN1, * IP1,DESCRN,ARN,IAN1,IAN2,INFON,IP2,LARN,LIAN1,
* LIAN2,AUX,LAUX,IERROR) * LIAN2,AUX,LAUX,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -59,8 +59,8 @@ C DESCRN Output matrix in JAD format
C ARN,IAN1 C ARN,IAN1
C IAN2,INFON, IP1, IP2 C IAN2,INFON, IP1, IP2
C C
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C C
C .. Scalar Arguments .. C .. Scalar Arguments ..

@ -31,8 +31,8 @@ C
SUBROUTINE ZGIND_TRI(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA, SUBROUTINE ZGIND_TRI(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA,
+ LARN,LKA,LJA,IPERM,WORK, LWORK, IERROR) + LARN,LKA,LJA,IPERM,WORK, LWORK, IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C ... Scalar arguments ... C ... Scalar arguments ...

@ -31,8 +31,8 @@ C
SUBROUTINE ZGINDEX(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA, INFON, SUBROUTINE ZGINDEX(M,N_BLOCKS,A,IA1,IA2,ARN,KA,IA,JA, INFON,
+ LARN,LKA,LJA,IPERM,WORK, LWORK, SIZE_REQ, IERROR) + LARN,LKA,LJA,IPERM,WORK, LWORK, SIZE_REQ, IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C ... Scalar arguments ... C ... Scalar arguments ...
INTEGER M, LWORK,N_BLOCKS,LARN,LKA,LJA, INTEGER M, LWORK,N_BLOCKS,LARN,LKA,LJA,

@ -190,8 +190,8 @@ C
C C
SUBROUTINE DCSMM(TRANS,M,N,K,ALPHA,PL,FIDA,DESCRA,A,IA1,IA2, SUBROUTINE DCSMM(TRANS,M,N,K,ALPHA,PL,FIDA,DESCRA,A,IA1,IA2,
& INFOA,PR,B,LDB,BETA,C,LDC,WORK,LWORK,IERROR) & INFOA,PR,B,LDB,BETA,C,LDC,WORK,LWORK,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C .. Scalar Arguments .. C .. Scalar Arguments ..
INTEGER M,N,K,LDB,LDC,LWORK, IERROR INTEGER M,N,K,LDB,LDC,LWORK, IERROR

@ -31,8 +31,8 @@ C
C ... Compute infinity norma for sparse matrix in CSR Format ... C ... Compute infinity norma for sparse matrix in CSR Format ...
DOUBLE PRECISION FUNCTION DJADNR(TRANS,M,N,NG,A,KA,JA,IA, DOUBLE PRECISION FUNCTION DJADNR(TRANS,M,N,NG,A,KA,JA,IA,
+ INFOA,IERROR) + INFOA,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C .. Scalar Arguments .. C .. Scalar Arguments ..
INTEGER M,N, IERROR, NG INTEGER M,N, IERROR, NG
CHARACTER TRANS CHARACTER TRANS

@ -31,8 +31,8 @@ C
C ... Compute infinity norma for sparse matrix in CSR Format ... C ... Compute infinity norma for sparse matrix in CSR Format ...
SUBROUTINE DJADRWS(TRANS,M,N,NG,A,KA,JA,IA, SUBROUTINE DJADRWS(TRANS,M,N,NG,A,KA,JA,IA,
+ INFOA,ROWSUM,IERROR) + INFOA,ROWSUM,IERROR)
use psb_const_mod
IMPLICIT NONE IMPLICIT NONE
INCLUDE 'psb_const.fh'
C .. Scalar Arguments .. C .. Scalar Arguments ..
INTEGER M,N, IERROR, NG INTEGER M,N, IERROR, NG
CHARACTER TRANS CHARACTER TRANS

@ -52,8 +52,8 @@ subroutine psb_cdasb(desc_a,info)
!....Locals.... !....Locals....
integer :: int_err(5), itemp(2) integer :: int_err(5), itemp(2)
integer,pointer :: ovrlap_index(:),halo_index(:),ovrlap_out(:),halo_out(:) integer,pointer :: ovrlap_index(:),halo_index(:)
integer :: i,err,np,me,& integer :: i,err,np,me,&
& lovrlap,lhalo,nhalo,novrlap,max_size,max_halo,n_col,ldesc_halo,& & lovrlap,lhalo,nhalo,novrlap,max_size,max_halo,n_col,ldesc_halo,&
& ldesc_ovrlap, dectype, err_act & ldesc_ovrlap, dectype, err_act
integer :: ictxt,n_row integer :: ictxt,n_row
@ -105,106 +105,28 @@ subroutine psb_cdasb(desc_a,info)
call psb_realloc(desc_a%matrix_data(psb_n_col_),desc_a%loc_to_glob,info) call psb_realloc(desc_a%matrix_data(psb_n_col_),desc_a%loc_to_glob,info)
! comm desc_size is size requested for temporary comm descriptors
! (expressed in No of dble element)
ldesc_halo = (((3*(n_col-n_row)+1)+1))
ovrlap_index => desc_a%ovrlap_index ovrlap_index => desc_a%ovrlap_index
nullify(desc_a%ovrlap_index,ovrlap_out) nullify(desc_a%ovrlap_index)
halo_index => desc_a%halo_index halo_index => desc_a%halo_index
nullify(desc_a%halo_index,halo_out) nullify(desc_a%halo_index)
lhalo = 1 call psi_cnv_dsc(halo_index,ovrlap_index,desc_a,info)
do while (halo_index(lhalo) /= -1) if (info /= 0) then
lhalo = lhalo + 1 call psb_errpush(4010,name,a_err='psi_bld_cdesc')
enddo
nhalo = (lhalo-1)/3
lovrlap=1
do while (ovrlap_index(lovrlap) /= -1)
lovrlap=lovrlap+1
enddo
novrlap = (lovrlap-1)/3
! Allocate final comm PSBLAS descriptors
! compute necessary dimension of halo index
max_halo = max(nhalo,1)
max_size = max(1,min(3*desc_a%matrix_data(psb_n_row_),novrlap*3))
itemp(1) = max_size
itemp(2) = max_halo
call psb_amx(ictxt, itemp(1:2))
max_size = itemp(1)
max_halo = itemp(2)
ldesc_halo = 3*max_halo+3*nhalo+1
! allocate HALO_INDEX field
call psb_realloc(ldesc_halo, halo_out, info)
! check on allocate
if (info /= no_err) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
! compute necessary dimension of ovrlap index
ldesc_ovrlap = 2*lovrlap+1
! allocate OVRLAP_INDEX field
call psb_realloc(ldesc_ovrlap, ovrlap_out, info)
! check on allocate
if (info /= no_err) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) write(0,*) 'psb_cdasb: converting indexes',&
& nhalo,lhalo,halo_index(lhalo)
!.... convert comunication stuctures....
! first the halo index
call psi_crea_index(desc_a,halo_index,&
& halo_out,.false.,info)
if(info.ne.0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
desc_a%halo_index => halo_out
! then the overlap index
call psi_crea_index(desc_a,ovrlap_index,&
& ovrlap_out,.true.,info)
if(info.ne.0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
desc_a%ovrlap_index => ovrlap_out
! next is the ovrlap_elem index
call psi_crea_ovr_elem(desc_a%ovrlap_index,desc_a%ovrlap_elem)
! finally bnd_elem
call psi_crea_bnd_elem(desc_a,info)
if(info.ne.0) then
call psb_errpush(4010,name,a_err='psi_crea_bnd_elem')
goto 9999 goto 9999
end if end if
! Ok, register into MATRIX_DATA & free temporary work areas ! Ok, register into MATRIX_DATA & free temporary work areas
desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_ desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_
deallocate(ovrlap_index, stat=info) deallocate(ovrlap_index, halo_index, stat=info)
deallocate(halo_index, stat=info)
if (info /= 0) then if (info /= 0) then
info =4000 info =4000
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
else else
info = 600 info = 600
call psb_errpush(info,name) call psb_errpush(info,name)

@ -54,7 +54,7 @@ subroutine psb_cdcpy(desc_in, desc_out, info)
integer, intent(out) :: info integer, intent(out) :: info
!locals !locals
integer :: np,me,ictxt, isz, err_act integer :: np,me,ictxt, isz, err_act, isz1, isz2
logical, parameter :: debug=.false.,debugprt=.false. logical, parameter :: debug=.false.,debugprt=.false.
character(len=20) :: name, char_err character(len=20) :: name, char_err
if (debug) write(0,*) me,'Entered CDCPY' if (debug) write(0,*) me,'Entered CDCPY'
@ -62,162 +62,37 @@ subroutine psb_cdcpy(desc_in, desc_out, info)
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_cdcpy' name = 'psb_cdcpy'
ictxt=desc_in%matrix_data(psb_ctxt_) ictxt=desc_in%matrix_data(psb_ctxt_)
! check on blacs grid ! check on blacs grid
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (debug) write(0,*) me,'Entered CDCPY' if (debug) write(0,*) me,'Entered CDCPY'
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
call psb_nullify_desc(desc_out) call psb_nullify_desc(desc_out)
if (associated(desc_in%matrix_data)) then call psb_safe_cpy(desc_in%matrix_data,desc_out%matrix_data,info)
isz = size(desc_in%matrix_data) if (info == 0) call psb_safe_cpy(desc_in%halo_index,desc_out%halo_index,info)
! allocate(desc_out%matrix_data(isz),stat=info) !!$ if (info == 0) call psb_safe_cpy(desc_in%halo_pt,desc_out%halo_pt,info)
call psb_realloc(isz,desc_out%matrix_data,info) if (info == 0) call psb_safe_cpy(desc_in%ovrlap_index,desc_out%ovrlap_index,info)
if(debug) write(0,*) 'cdcpy: m_data',isz,':',desc_in%matrix_data(:) !!$ if (info == 0) call psb_safe_cpy(desc_in%ovrlap_pt,desc_out%ovrlap_pt,info)
if (info /= 0) then if (info == 0) call psb_safe_cpy(desc_in%bnd_elem,desc_out%bnd_elem,info)
info=4010 if (info == 0) call psb_safe_cpy(desc_in%ovrlap_elem,desc_out%ovrlap_elem,info)
char_err='psb_realloc' if (info == 0) call psb_safe_cpy(desc_in%loc_to_glob,desc_out%loc_to_glob,info)
call psb_errpush(info,name,a_err=char_err) if (info == 0) call psb_safe_cpy(desc_in%glob_to_loc,desc_out%glob_to_loc,info)
goto 9999 if (info == 0) call psb_safe_cpy(desc_in%lprm,desc_out%lprm,info)
else if (info == 0) call psb_safe_cpy(desc_in%idx_space,desc_out%idx_space,info)
desc_out%matrix_data(:) = desc_in%matrix_data(:)
endif
endif
if (debug) write(0,*) me,'Done matrix_data '
if (associated(desc_in%halo_index)) then
isz = size(desc_in%halo_index)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%halo_index,info)
if(debugprt) write(0,*) 'cdcpy: h_idx',isz,':',desc_in%halo_index(:)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%halo_index(:) = desc_in%halo_index(:)
endif
endif
if (debug) write(0,*) me,'Done halo_index'
if (associated(desc_in%bnd_elem)) then if (info /= 0) then
isz = size(desc_in%bnd_elem) info = 4010
! allocate(desc_out%matrix_data(isz),stat=info) call psb_errpush(info,name)
call psb_realloc(isz,desc_out%bnd_elem,info) goto 9999
if(debugprt) write(0,*) 'cdcpy: bnd_elem',isz,':',desc_in%bnd_elem(:)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%bnd_elem(:) = desc_in%bnd_elem(:)
endif
endif endif
if (debug) write(0,*) me,'Done bnd_elem'
if (associated(desc_in%ovrlap_elem)) then
isz = size(desc_in%ovrlap_elem)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%ovrlap_elem,info)
if(debugprt) write(0,*) 'cdcpy: ovrlap_elem',isz,':',desc_in%ovrlap_elem(:)
if (info /= 0) then
info=4010
char_err='psrealloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%ovrlap_elem(:) = desc_in%ovrlap_elem(:)
endif
endif
if (debug) write(0,*) me,'Done ovrlap_elem'
if (associated(desc_in%ovrlap_index)) then
isz = size(desc_in%ovrlap_index)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%ovrlap_index,info)
if(debugprt) write(0,*) 'cdcpy: ovrlap_index',isz,':',desc_in%ovrlap_index(:)
if (info /= 0) then
info=4010
char_err='psrealloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%ovrlap_index(:) = desc_in%ovrlap_index(:)
endif
endif
if (debug) write(0,*) me,'Done ovrlap_index'
if (associated(desc_in%loc_to_glob)) then
isz = size(desc_in%loc_to_glob)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%loc_to_glob,info)
if(debugprt) write(0,*) 'cdcpy: loc_to_glob',isz,':',desc_in%loc_to_glob(:)
if (info /= 0) then
info=4010
char_err='psrealloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%loc_to_glob(:) = desc_in%loc_to_glob(:)
endif
endif
if (debug) write(0,*) me,'Done loc_to_glob'
if (associated(desc_in%glob_to_loc)) then
isz = size(desc_in%glob_to_loc)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%glob_to_loc,info)
if(debugprt) write(0,*) 'cdcpy: glob_to_loc',isz,':',desc_in%glob_to_loc(:)
if (info /= 0) then
info=4010
char_err='psrealloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%glob_to_loc(:) = desc_in%glob_to_loc(:)
endif
endif
if (debug) write(0,*) me,'Done glob_to_loc'
if (associated(desc_in%lprm)) then
isz = size(desc_in%lprm)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%lprm,info)
if(debugprt) write(0,*) 'cdcpy: lprm',isz,':',desc_in%lprm(:)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%lprm(:) = desc_in%lprm(:)
endif
endif
if (debug) write(0,*) me,'Done lprm'
if (associated(desc_in%idx_space)) then
isz = size(desc_in%idx_space)
! allocate(desc_out%matrix_data(isz),stat=info)
call psb_realloc(isz,desc_out%idx_space,info)
if(debugprt) write(0,*) 'cdcpy: idx_space',isz,':',desc_in%idx_space(:)
if (info /= 0) then
info=4010
char_err='psb_realloc'
call psb_errpush(info,name,a_err=char_err)
goto 9999
else
desc_out%idx_space(:) = desc_in%idx_space(:)
endif
endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -107,6 +107,7 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
use psb_serial_mod use psb_serial_mod
use psb_const_mod use psb_const_mod
use psb_error_mod use psb_error_mod
use psi_mod
use psb_penv_mod use psb_penv_mod
implicit None implicit None
!....Parameters... !....Parameters...
@ -117,7 +118,7 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
!locals !locals
Integer :: i,j,np,me,err,n,itmpov, k,& Integer :: i,j,np,me,err,n,itmpov, k,&
& l_ov_ix,l_ov_el,idx, err_act,m, ip & l_ov_ix,l_ov_el,idx, err_act,m, ip
Integer :: INT_ERR(5) Integer :: INT_ERR(5), thalo(1), tovr(1)
integer, allocatable :: nlv(:) integer, allocatable :: nlv(:)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name character(len=20) :: name
@ -164,10 +165,12 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
!count local rows number !count local rows number
! allocate work vector ! allocate work vector
!!$ allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),&
!!$ & desc_a%loc_to_glob(nloc),desc_a%lprm(1),&
!!$ & desc_a%ovrlap_index(1),desc_a%ovrlap_elem(1),&
!!$ & desc_a%halo_index(1),desc_a%bnd_elem(1),stat=info)
allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),& allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),&
& desc_a%loc_to_glob(nloc),desc_a%lprm(1),& & desc_a%loc_to_glob(m),desc_a%lprm(1),stat=info)
& desc_a%ovrlap_index(1),desc_a%ovrlap_elem(1),&
& desc_a%halo_index(1),desc_a%bnd_elem(1),stat=info)
if (info /= 0) then if (info /= 0) then
info=2025 info=2025
int_err(1)=m int_err(1)=m
@ -175,7 +178,13 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
goto 9999 goto 9999
endif endif
desc_a%matrix_data(psb_m_) = m
desc_a%matrix_data(psb_n_) = m
desc_a%matrix_data(psb_n_row_) = nloc
desc_a%matrix_data(psb_n_col_) = nloc
desc_a%matrix_data(psb_dec_type_) = psb_desc_bld_
desc_a%matrix_data(psb_ctxt_) = ictxt
call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_))
j = 1 j = 1
do ip=0, np-1 do ip=0, np-1
@ -194,21 +203,18 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
enddo enddo
tovr = -1
thalo = -1
desc_a%lprm(:) = 0 desc_a%lprm(:) = 0
desc_a%halo_index(:) = -1
desc_a%bnd_elem(:) = -1
desc_a%ovrlap_index(:) = -1
desc_a%ovrlap_elem(:) = -1
call psi_cnv_dsc(thalo,tovr,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psi_bld_cdesc')
goto 9999
end if
desc_a%matrix_data(psb_m_) = m
desc_a%matrix_data(psb_n_) = m
desc_a%matrix_data(psb_n_row_) = nloc
desc_a%matrix_data(psb_n_col_) = nloc
desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_ desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_
desc_a%matrix_data(psb_ctxt_) = ictxt
call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_))
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -159,6 +159,34 @@ subroutine psb_cdfree(desc_a,info)
goto 9999 goto 9999
end if end if
if (associated(desc_a%idx_space)) then
deallocate(desc_a%idx_space,stat=info)
if (info /= 0) then
info=2056
call psb_errpush(info,name)
goto 9999
end if
end if
!!$ if (associated(desc_a%halo_pt)) then
!!$ deallocate(desc_a%halo_pt,stat=info)
!!$ if (info /= 0) then
!!$ info=2056
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!!$ end if
!!$
!!$ if (associated(desc_a%ovrlap_pt)) then
!!$ deallocate(desc_a%ovrlap_pt,stat=info)
!!$ if (info /= 0) then
!!$ info=2056
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!!$ end if
call psb_nullify_desc(desc_a) call psb_nullify_desc(desc_a)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

@ -62,7 +62,7 @@ subroutine psb_cdprt(iout,desc_p,glob,short)
endif endif
if (.not.lglob) then if (.not.lglob) then
write(iout,*) 'Precond. descriptor:',desc_p%matrix_data(1:10) write(iout,*) 'Communication descriptor:',desc_p%matrix_data(1:10)
m=desc_p%matrix_data(psb_m_) m=desc_p%matrix_data(psb_m_)
n_row=desc_p%matrix_data(psb_n_row_) n_row=desc_p%matrix_data(psb_n_row_)
n_col=desc_p%matrix_data(psb_n_col_) n_col=desc_p%matrix_data(psb_n_col_)
@ -121,7 +121,7 @@ subroutine psb_cdprt(iout,desc_p,glob,short)
else if (lglob) then else if (lglob) then
write(iout,*) 'Precond. descriptor:',desc_p%matrix_data(1:10) write(iout,*) 'Communication descriptor:',desc_p%matrix_data(1:10)
m=desc_p%matrix_data(psb_m_) m=desc_p%matrix_data(psb_m_)
n_row=desc_p%matrix_data(psb_n_row_) n_row=desc_p%matrix_data(psb_n_row_)
n_col=desc_p%matrix_data(psb_n_col_) n_col=desc_p%matrix_data(psb_n_col_)

@ -107,6 +107,7 @@ subroutine psb_cdrep(m, ictxt, desc_a, info)
use psb_serial_mod use psb_serial_mod
use psb_const_mod use psb_const_mod
use psb_error_mod use psb_error_mod
use psi_mod
use psb_penv_mod use psb_penv_mod
implicit None implicit None
!....Parameters... !....Parameters...
@ -116,7 +117,7 @@ subroutine psb_cdrep(m, ictxt, desc_a, info)
!locals !locals
Integer :: i,np,me,err,n,err_act Integer :: i,np,me,err,n,err_act
integer :: int_err(5),exch(2) integer :: int_err(5),exch(2), thalo(1), tovr(1)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name character(len=20) :: name
@ -174,9 +175,9 @@ subroutine psb_cdrep(m, ictxt, desc_a, info)
!count local rows number !count local rows number
! allocate work vector ! allocate work vector
allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),& allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),&
& desc_a%loc_to_glob(m),desc_a%lprm(1),& & desc_a%loc_to_glob(m),desc_a%lprm(1),stat=info)
& desc_a%ovrlap_index(1),desc_a%ovrlap_elem(1),& !!$ & desc_a%ovrlap_index(1),desc_a%ovrlap_elem(1),&
& desc_a%halo_index(1),desc_a%bnd_elem(1),stat=info) !!$ & desc_a%halo_index(1),desc_a%bnd_elem(1),stat=info)
if (info /= 0) then if (info /= 0) then
info=2025 info=2025
int_err(1)=m int_err(1)=m
@ -184,25 +185,32 @@ subroutine psb_cdrep(m, ictxt, desc_a, info)
goto 9999 goto 9999
endif endif
desc_a%matrix_data(psb_m_) = m
desc_a%matrix_data(psb_n_) = n
desc_a%matrix_data(psb_n_row_) = m
desc_a%matrix_data(psb_n_col_) = n
desc_a%matrix_data(psb_dec_type_) = psb_desc_bld_
desc_a%matrix_data(psb_ctxt_) = ictxt
call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_))
do i=1,m do i=1,m
desc_a%glob_to_loc(i) = i desc_a%glob_to_loc(i) = i
desc_a%loc_to_glob(i) = i desc_a%loc_to_glob(i) = i
enddo enddo
desc_a%lprm(:) = 0
desc_a%halo_index(:) = -1
desc_a%bnd_elem(:) = -1
desc_a%ovrlap_index(:) = -1
desc_a%ovrlap_elem(:) = -1
desc_a%matrix_data(psb_m_) = m tovr = -1
desc_a%matrix_data(psb_n_) = n thalo = -1
desc_a%matrix_data(psb_n_row_) = m desc_a%lprm(:) = 0
desc_a%matrix_data(psb_n_col_) = n
call psi_cnv_dsc(thalo,tovr,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psi_bld_cdesc')
goto 9999
end if
desc_a%matrix_data(psb_dec_type_) = psb_desc_repl_ desc_a%matrix_data(psb_dec_type_) = psb_desc_repl_
desc_a%matrix_data(psb_ctxt_) = ictxt
call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_))
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -170,10 +170,9 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
index_dim = size(desc_a%halo_index) index_dim = size(desc_a%halo_index)
elem_dim = size(desc_a%halo_index) elem_dim = size(desc_a%halo_index)
allocate(desc_ov%ovrlap_index(novr*(Max(2*index_dim,1)+1)),& nullify(desc_ov%ovrlap_index,desc_ov%halo_index,desc_ov%ovrlap_elem)
& desc_ov%ovrlap_elem(novr*(Max(elem_dim,1)+3)),& allocate(desc_ov%ovrlap_elem(novr*(Max(elem_dim,1)+3)),&
& desc_ov%matrix_data(psb_mdata_size_),& & desc_ov%matrix_data(psb_mdata_size_),STAT=INFO)
& desc_ov%halo_index(novr*(Size(desc_a%halo_index)+3)),STAT=INFO)
if (info /= 0) then if (info /= 0) then
info=4000 info=4000
call psb_errpush(info,name) call psb_errpush(info,name)
@ -183,10 +182,7 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
l_tmp_ovr_idx=novr*(3*Max(2*index_dim,1)+1) l_tmp_ovr_idx=novr*(3*Max(2*index_dim,1)+1)
l_tmp_halo=novr*(3*Size(desc_a%halo_index)) l_tmp_halo=novr*(3*Size(desc_a%halo_index))
desc_ov%ovrlap_index(:) = -1 desc_ov%matrix_data(:) = desc_a%matrix_data(:)
desc_ov%ovrlap_elem(:) = -1
desc_ov%halo_index(:) = -1
desc_ov%matrix_data(1:10) = desc_a%matrix_data(1:10)
desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_
Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),& Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),&
@ -203,11 +199,13 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
Write(0,*)'Start cdovrbld',me,lworks,lworkr Write(0,*)'Start cdovrbld',me,lworks,lworkr
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif endif
! !
! The real work goes on in here.... ! The real work goes on in here....
! !
Call psb_cdovrbld(novr,desc_ov,desc_a,a,& Call psb_cdovrbld(novr,desc_ov,desc_a,a,&
& l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psb_cdovrbld' ch_err='psb_cdovrbld'

@ -38,14 +38,16 @@
! !
! Parameters: ! Parameters:
! n_ovr - integer. The number of overlap levels ! n_ovr - integer. The number of overlap levels
! desc_p - type(<psb_desc_type>). The communication descriptor for the preconditioner. ! desc_p - type(<psb_desc_type>). The communication descriptor
! for the preconditioner.
! desc_a - type(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! a - type(<psb_dspmat_type). The matrix upon which the preconditioner will be built. ! a - type(<psb_dspmat_type). The matrix upon which the preconditioner
! will be built.
! l_tmp_halo - integer. Input estimate for allocation sizes. ! l_tmp_halo - integer. Input estimate for allocation sizes.
! l_tmp_ovr_idx - integer. Input estimate for allocation sizes. ! l_tmp_ovr_idx - integer. Input estimate for allocation sizes.
! lworkr - integer. Input estimate for allocation sizes. ! lworkr - integer. Input estimate for allocation sizes.
! lworks - integer. Input estimate for allocation sizes. ! lworks - integer. Input estimate for allocation sizes.
! info - integer. Eventually returns an error code ! info - integer. Possibly return an error code
Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
& l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info)
use psb_descriptor_type use psb_descriptor_type
@ -76,12 +78,11 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,& Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,&
& n_elem_send,tot_recv,tot_elem,n_col,m,ictxt,np,me,dl_lda,lwork,& & n_elem_send,tot_recv,tot_elem,n_col,m,ictxt,np,me,dl_lda,lwork,&
& counter_t,n_elem,i_ovr,jj,i,proc_id,isz, mglob, glx,n_row, & & counter_t,n_elem,i_ovr,jj,i,proc_id,isz, mglob, glx,n_row, &
& idxr, idxs, lx, iszr, err_act, icomm & idxr, idxs, lx, iszr, err_act, icomm, nxch, nsnd, nrcv
Integer,Pointer :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),& Integer,Pointer :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),&
& t_halo_out(:),work(:),dep_list(:),temp(:) & t_halo_out(:),work(:),dep_list(:),temp(:)
Integer,Pointer :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) Integer,Pointer :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:)
integer :: pairtree(2)
Logical,Parameter :: debug=.false. Logical,Parameter :: debug=.false.
real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7, tl, tch real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7, tl, tch
@ -117,8 +118,8 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
dl_lda=np*5 dl_lda=np*5
lwork=5*(5*np+2)*np+10 lwork=5*(5*np+2)*np+10
Allocate(works(lworks),workr(lworkr),t_halo_in(3*Size(desc_p%halo_index)),& Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),&
& t_halo_out(Size(desc_p%halo_index)), work(lwork),& & t_halo_out(l_tmp_halo), work(lwork),&
& length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info) & length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -152,14 +153,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
counter_o = 1 counter_o = 1
! See comment in main loop below. ! See comment in main loop below.
call InitPairSearchTree(pairtree,info)
if (info /= 0) then
info=4010
ch_err='InitPairSearhTree'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) me,'Done InitPairSearchTree',info
! !
! A picture is in order to understand what goes on here. ! A picture is in order to understand what goes on here.
@ -216,16 +209,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
tot_recv=tot_recv+n_elem_recv tot_recv=tot_recv+n_elem_recv
if (debug) write(0,*) me,' CDOVRBLD tot_recv:',proc,n_elem_recv,tot_recv if (debug) write(0,*) me,' CDOVRBLD tot_recv:',proc,n_elem_recv,tot_recv
! !
! While running through the column indices exchanged with other procs
! we have to keep track of which elements actually are overlap and halo
! ones to record them in overlap_elem. We do this by maintaining
! an AVL balanced search tree: at each point counter_e is the next
! free index element. The search routine for gidx will return
! glx if gidx was already assigned a local index (glx<counter_e)
! but if gidx was a new index for this process, then it creates
! a new pair (gidx,counter_e), and glx==counter_e. In this case we
! need to record this for the overlap exchange. In the first iteration
! this gets filled with the first halo indices.
! !
! The format of the halo vector exists in two forms: 1. Temporary ! The format of the halo vector exists in two forms: 1. Temporary
! 2. Assembled. In this loop we are using the (assembled) halo_in and ! 2. Assembled. In this loop we are using the (assembled) halo_in and
@ -234,13 +217,13 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
! everything for the next iteration. ! everything for the next iteration.
! !
!!$ if (me==0) Write(0,*)'Loop ',size(halo), counter, n_elem_recv,n_elem_send
t3 = mpi_wtime() t3 = mpi_wtime()
! !
! add recv elements in halo_index into ovrlap_index ! add recv elements in halo_index into ovrlap_index
! !
Do j=0,n_elem_recv-1 Do j=0,n_elem_recv-1
If((counter+psb_elem_recv_+j)>Size(halo)) then
If ((counter+psb_elem_recv_+j)>Size(halo)) then
info=-2 info=-2
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -292,34 +275,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
counter_h=counter_h+3 counter_h=counter_h+3
call SearchInsKeyVal(pairtree,gidx,counter_e,glx,info)
!!$ if (debug) write(0,*) 'From searchInsKey ',gidx,glx,counter_e,info
if (info>=0) then
If (glx < counter_e) Then
desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)= &
& desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)+1
Else
If((counter_e+2) > Size(desc_p%ovrlap_elem)) Then
isz = max((3*Size(desc_p%ovrlap_elem))/2,(counter_e+3))
if (debug) write(0,*) me,'Realloc ovr_El',isz
call psb_realloc(isz,desc_p%ovrlap_elem,info)
if (info /= 0) then
info=4010
ch_err='psrealloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
!!$ if (debug) write(0,*) 'Adding into ovrlap ',gidx,glx,counter_e,info
desc_p%ovrlap_elem(counter_e)=gidx
desc_p%ovrlap_elem(counter_e+psb_n_dom_ovr_)=2
desc_p%ovrlap_elem(counter_e+2)=-1
counter_e = counter_e + 2
End If
else
write(0,*) me, 'Cdovrbld From SearchInsKeyVal: ',info
endif
Enddo Enddo
if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10)
counter = counter+n_elem_recv counter = counter+n_elem_recv
@ -350,34 +305,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
tmp_ovr_idx(counter_o+3)=-1 tmp_ovr_idx(counter_o+3)=-1
counter_o=counter_o+3 counter_o=counter_o+3
call SearchInsKeyVal(pairtree,gidx,counter_e,glx,info)
!!$ if (debug) write(0,*) 'From searchInsKey ',gidx,glx,counter_e,info
if (info>=0) then
If (glx < counter_e) Then
desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)= &
& desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)+1
Else
If((counter_e+2) > Size(desc_p%ovrlap_elem)) Then
isz = max((3*Size(desc_p%ovrlap_elem))/2,(counter_e+3))
if (debug) write(0,*) me,'Realloc ovr_el',isz
call psb_realloc(isz,desc_p%ovrlap_elem,info)
if (info /= 0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
!!$ if (debug) write(0,*) 'Adding into ovrlap ',gidx,glx,counter_e,info
desc_p%ovrlap_elem(counter_e)=gidx
desc_p%ovrlap_elem(counter_e+psb_n_dom_ovr_)=2
desc_p%ovrlap_elem(counter_e+2)=-1
counter_e = counter_e + 2
End If
else
write(0,*) me,'Cdovrbld From SearchInsKeyVal: ',info
endif
! !
! Prepare to exchange the halo rows with the other proc. ! Prepare to exchange the halo rows with the other proc.
! !
@ -402,6 +329,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
end if end if
lworks = isz lworks = isz
End If End If
If((n_elem) > size(blk%ia2)) Then If((n_elem) > size(blk%ia2)) Then
isz = max((3*size(blk%ia2))/2,(n_elem)) isz = max((3*size(blk%ia2))/2,(n_elem))
if (debug) write(0,*) me,'Realloc blk',isz if (debug) write(0,*) me,'Realloc blk',isz
@ -434,18 +362,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
if (i_ovr < n_ovr) then if (i_ovr < n_ovr) then
if (tot_elem > 1) then if (tot_elem > 1) then
!!$ write(0,*) me,'Realloc temp',tot_elem+2 call imsr(tot_elem,works(idxs+1))
if (tot_elem+2 > size(temp)) then
!!$ write(0,*) me,'Realloc temp',tot_elem+2
deallocate(temp)
allocate(temp(tot_elem+2),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
Call mrgsrt(tot_elem,works(idxs+1),temp,info)
If (info == 0) Call ireordv1(tot_elem,works(idxs+1),temp)
lx = works(idxs+1) lx = works(idxs+1)
i = 1 i = 1
j = 1 j = 1
@ -546,7 +463,18 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
End If End If
desc_p%glob_to_loc(idx)=n_col desc_p%glob_to_loc(idx)=n_col
desc_p%loc_to_glob(n_col)=idx desc_p%loc_to_glob(n_col)=idx
If((counter_t+3) > Size(t_halo_in))Write(0,*)'bingo' If ((counter_t+3) > Size(t_halo_in)) then
isz = max((3*Size(t_halo_in))/2,(counter_t+3+1000))
if (debug) write(0,*) me,'Realloc ovr_el',isz
call psb_realloc(isz,t_halo_in,info)
if (info /= 0) then
info=4010
ch_err='psrealloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end If
t_halo_in(counter_t)=proc_id t_halo_in(counter_t)=proc_id
t_halo_in(counter_t+1)=1 t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col t_halo_in(counter_t+2)=n_col
@ -591,7 +519,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Calling Crea_Halo' if (debug) write(0,*) me,'Calling Crea_Halo'
call psi_crea_index(desc_p,t_halo_in,t_halo_out,.false.,info) call psi_crea_index(desc_p,t_halo_in,t_halo_out,.false.,nxch,nsnd,nrcv,info)
if (debug) then if (debug) then
write(0,*) me,'Done Crea_Index' write(0,*) me,'Done Crea_Index'
@ -611,7 +539,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
tch = tch +(t3-t2) tch = tch +(t3-t2)
End Do End Do
t1 = mpi_wtime() t1 = mpi_wtime()
call FreePairSearchTree(pairtree)
desc_p%matrix_data(psb_m_)=desc_a%matrix_data(psb_m_) desc_p%matrix_data(psb_m_)=desc_a%matrix_data(psb_m_)
desc_p%matrix_data(psb_n_)=desc_a%matrix_data(psb_n_) desc_p%matrix_data(psb_n_)=desc_a%matrix_data(psb_n_)
@ -627,35 +554,12 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
! !
if (debug) then if (debug) then
write(0,*) 'psb_cdasb: converting indexes' write(0,*) 'psb_cdovrbld: converting indexes'
call psb_barrier(ictxt) call psb_barrier(ictxt)
end if end if
!.... convert comunication stuctures.... !.... convert comunication stuctures....
! first the halo index
call psi_crea_index(desc_p,tmp_halo,&
& desc_p%halo_index,.false.,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
! then the overlap index
call psi_crea_index(desc_p,tmp_ovr_idx,&
& desc_p%ovrlap_index,.true.,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
! next is the ovrlap_elem index call psi_cnv_dsc(tmp_halo,tmp_ovr_idx,desc_p,info)
call psi_crea_ovr_elem(desc_p%ovrlap_index,desc_p%ovrlap_elem)
! finally bnd_elem
call psi_crea_bnd_elem(desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_bnd_elem')
goto 9999
end if
! Ok, register into MATRIX_DATA & free temporary work areas ! Ok, register into MATRIX_DATA & free temporary work areas
desc_p%matrix_data(psb_dec_type_) = psb_desc_asb_ desc_p%matrix_data(psb_dec_type_) = psb_desc_asb_

@ -54,10 +54,9 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
use psb_tools_mod, only : psb_glob_to_loc, psb_loc_to_glob use psb_tools_mod, only : psb_glob_to_loc, psb_loc_to_glob
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use mpi
Implicit None Implicit None
include 'mpif.h'
Type(psb_dspmat_type),Intent(in) :: a Type(psb_dspmat_type),Intent(in) :: a
Type(psb_dspmat_type),Intent(inout) :: blk Type(psb_dspmat_type),Intent(inout) :: blk
Type(psb_desc_type),Intent(in) :: desc_a Type(psb_desc_type),Intent(in) :: desc_a

@ -91,7 +91,8 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
! .. Local Scalars .. ! .. Local Scalars ..
Integer :: np, me,m,nnzero, ictxt, lovr, lworks,lworkr, n_col, int_err(5),& Integer :: i, j, k, np, me,m,nnzero,&
& ictxt, lovr, lworks,lworkr, n_col, int_err(5),&
& index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo & index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo
Logical, parameter :: debug=.false. Logical, parameter :: debug=.false.
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -144,7 +145,6 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
t1 = mpi_wtime() t1 = mpi_wtime()
! !
! Ok, since we are only estimating, do it as follows: ! Ok, since we are only estimating, do it as follows:
! LOVR= (NNZ/NROW)*N_HALO*N_OVR This assumes that the local average ! LOVR= (NNZ/NROW)*N_HALO*N_OVR This assumes that the local average
@ -170,10 +170,9 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
index_dim = size(desc_a%halo_index) index_dim = size(desc_a%halo_index)
elem_dim = size(desc_a%halo_index) elem_dim = size(desc_a%halo_index)
allocate(desc_ov%ovrlap_index(novr*(Max(2*index_dim,1)+1)),& nullify(desc_ov%ovrlap_index,desc_ov%halo_index,desc_ov%ovrlap_elem)
& desc_ov%ovrlap_elem(novr*(Max(elem_dim,1)+3)),& allocate(desc_ov%ovrlap_elem(novr*(Max(elem_dim,1)+3)),&
& desc_ov%matrix_data(psb_mdata_size_),& & desc_ov%matrix_data(psb_mdata_size_),STAT=INFO)
& desc_ov%halo_index(novr*(Size(desc_a%halo_index)+3)),STAT=INFO)
if (info /= 0) then if (info /= 0) then
info=4000 info=4000
call psb_errpush(info,name) call psb_errpush(info,name)
@ -183,10 +182,7 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
l_tmp_ovr_idx=novr*(3*Max(2*index_dim,1)+1) l_tmp_ovr_idx=novr*(3*Max(2*index_dim,1)+1)
l_tmp_halo=novr*(3*Size(desc_a%halo_index)) l_tmp_halo=novr*(3*Size(desc_a%halo_index))
desc_ov%ovrlap_index(:) = -1 desc_ov%matrix_data(:) = desc_a%matrix_data(:)
desc_ov%ovrlap_elem(:) = -1
desc_ov%halo_index(:) = -1
desc_ov%matrix_data(1:10) = desc_a%matrix_data(1:10)
desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_
Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),& Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),&

@ -37,15 +37,17 @@
! Note that n_ovr > 0 thanks to the caller routine. ! Note that n_ovr > 0 thanks to the caller routine.
! !
! Parameters: ! Parameters:
! n_ovr - integer. The number of overlap levels ! n_ovr - integer. The number of overlap levels
! desc_p - type(<psb_desc_type>). The communication descriptor for the preconditioner. ! desc_p - type(<psb_desc_type>). The communication descriptor
! for the preconditioner.
! desc_a - type(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! a - type(<psb_dspmat_type). The matrix upon which the preconditioner will be built. ! a - type(<psb_dspmat_type). The matrix upon which the preconditioner
! will be built.
! l_tmp_halo - integer. Input estimate for allocation sizes. ! l_tmp_halo - integer. Input estimate for allocation sizes.
! l_tmp_ovr_idx - integer. Input estimate for allocation sizes. ! l_tmp_ovr_idx - integer. Input estimate for allocation sizes.
! lworkr - integer. Input estimate for allocation sizes. ! lworkr - integer. Input estimate for allocation sizes.
! lworks - integer. Input estimate for allocation sizes. ! lworks - integer. Input estimate for allocation sizes.
! info - integer. Eventually returns an error code ! info - integer. Possibly returns an error code
Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
& l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info)
use psb_descriptor_type use psb_descriptor_type
@ -76,12 +78,11 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,& Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,&
& n_elem_send,tot_recv,tot_elem,n_col,m,ictxt,np,me,dl_lda,lwork,& & n_elem_send,tot_recv,tot_elem,n_col,m,ictxt,np,me,dl_lda,lwork,&
& counter_t,n_elem,i_ovr,jj,i,proc_id,isz, mglob, glx,n_row, & & counter_t,n_elem,i_ovr,jj,i,proc_id,isz, mglob, glx,n_row, &
& idxr, idxs, lx, iszr, err_act, icomm & idxr, idxs, lx, iszr, err_act, icomm,nxch,nsnd,nrcv
Integer,Pointer :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),& Integer,Pointer :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),&
& t_halo_out(:),work(:),dep_list(:),temp(:) & t_halo_out(:),work(:),dep_list(:),temp(:)
Integer,Pointer :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) Integer,Pointer :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:)
integer :: pairtree(2)
Logical,Parameter :: debug=.false. Logical,Parameter :: debug=.false.
real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7, tl, tch real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7, tl, tch
@ -117,8 +118,8 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
dl_lda=np*5 dl_lda=np*5
lwork=5*(5*np+2)*np+10 lwork=5*(5*np+2)*np+10
Allocate(works(lworks),workr(lworkr),t_halo_in(3*Size(desc_p%halo_index)),& Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),&
& t_halo_out(Size(desc_p%halo_index)), work(lwork),& & t_halo_out(l_tmp_halo), work(lwork),&
& length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info) & length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
@ -152,14 +153,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
counter_o = 1 counter_o = 1
! See comment in main loop below. ! See comment in main loop below.
call InitPairSearchTree(pairtree,info)
if (info /= 0) then
info=4010
ch_err='InitPairSearhTree'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) me,'Done InitPairSearchTree',info
! !
! A picture is in order to understand what goes on here. ! A picture is in order to understand what goes on here.
@ -184,7 +177,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
Do i_ovr=1,n_ovr Do i_ovr=1,n_ovr
if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',n_ovr if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',n_ovr
!!$ t_halo_in(:) = -1
! !
! At this point, halo contains a valid halo corresponding to the ! At this point, halo contains a valid halo corresponding to the
@ -216,16 +208,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
tot_recv=tot_recv+n_elem_recv tot_recv=tot_recv+n_elem_recv
if (debug) write(0,*) me,' CDOVRBLD tot_recv:',proc,n_elem_recv,tot_recv if (debug) write(0,*) me,' CDOVRBLD tot_recv:',proc,n_elem_recv,tot_recv
! !
! While running through the column indices exchanged with other procs
! we have to keep track of which elements actually are overlap and halo
! ones to record them in overlap_elem. We do this by maintaining
! an AVL balanced search tree: at each point counter_e is the next
! free index element. The search routine for gidx will return
! glx if gidx was already assigned a local index (glx<counter_e)
! but if gidx was a new index for this process, then it creates
! a new pair (gidx,counter_e), and glx==counter_e. In this case we
! need to record this for the overlap exchange. In the first iteration
! this gets filled with the first halo indices.
! !
! The format of the halo vector exists in two forms: 1. Temporary ! The format of the halo vector exists in two forms: 1. Temporary
! 2. Assembled. In this loop we are using the (assembled) halo_in and ! 2. Assembled. In this loop we are using the (assembled) halo_in and
@ -234,13 +216,13 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
! everything for the next iteration. ! everything for the next iteration.
! !
!!$ if (me==0) Write(0,*)'Loop ',size(halo), counter, n_elem_recv,n_elem_send
t3 = mpi_wtime() t3 = mpi_wtime()
! !
! add recv elements in halo_index into ovrlap_index ! add recv elements in halo_index into ovrlap_index
! !
Do j=0,n_elem_recv-1 Do j=0,n_elem_recv-1
If((counter+psb_elem_recv_+j)>Size(halo)) then
If ((counter+psb_elem_recv_+j)>Size(halo)) then
info=-2 info=-2
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -292,34 +274,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
counter_h=counter_h+3 counter_h=counter_h+3
call SearchInsKeyVal(pairtree,gidx,counter_e,glx,info)
!!$ if (debug) write(0,*) 'From searchInsKey ',gidx,glx,counter_e,info
if (info>=0) then
If (glx < counter_e) Then
desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)= &
& desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)+1
Else
If((counter_e+2) > Size(desc_p%ovrlap_elem)) Then
isz = max((3*Size(desc_p%ovrlap_elem))/2,(counter_e+3))
if (debug) write(0,*) me,'Realloc ovr_El',isz
call psb_realloc(isz,desc_p%ovrlap_elem,info)
if (info /= 0) then
info=4010
ch_err='psrealloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
!!$ if (debug) write(0,*) 'Adding into ovrlap ',gidx,glx,counter_e,info
desc_p%ovrlap_elem(counter_e)=gidx
desc_p%ovrlap_elem(counter_e+psb_n_dom_ovr_)=2
desc_p%ovrlap_elem(counter_e+2)=-1
counter_e = counter_e + 2
End If
else
write(0,*) me, 'Cdovrbld From SearchInsKeyVal: ',info
endif
Enddo Enddo
if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10)
counter = counter+n_elem_recv counter = counter+n_elem_recv
@ -350,34 +304,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
tmp_ovr_idx(counter_o+3)=-1 tmp_ovr_idx(counter_o+3)=-1
counter_o=counter_o+3 counter_o=counter_o+3
call SearchInsKeyVal(pairtree,gidx,counter_e,glx,info)
!!$ if (debug) write(0,*) 'From searchInsKey ',gidx,glx,counter_e,info
if (info>=0) then
If (glx < counter_e) Then
desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)= &
& desc_p%ovrlap_elem(glx+psb_n_dom_ovr_)+1
Else
If((counter_e+2) > Size(desc_p%ovrlap_elem)) Then
isz = max((3*Size(desc_p%ovrlap_elem))/2,(counter_e+3))
if (debug) write(0,*) me,'Realloc ovr_el',isz
call psb_realloc(isz,desc_p%ovrlap_elem,info)
if (info /= 0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
!!$ if (debug) write(0,*) 'Adding into ovrlap ',gidx,glx,counter_e,info
desc_p%ovrlap_elem(counter_e)=gidx
desc_p%ovrlap_elem(counter_e+psb_n_dom_ovr_)=2
desc_p%ovrlap_elem(counter_e+2)=-1
counter_e = counter_e + 2
End If
else
write(0,*) me,'Cdovrbld From SearchInsKeyVal: ',info
endif
! !
! Prepare to exchange the halo rows with the other proc. ! Prepare to exchange the halo rows with the other proc.
! !
@ -402,6 +328,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
end if end if
lworks = isz lworks = isz
End If End If
If((n_elem) > size(blk%ia2)) Then If((n_elem) > size(blk%ia2)) Then
isz = max((3*size(blk%ia2))/2,(n_elem)) isz = max((3*size(blk%ia2))/2,(n_elem))
if (debug) write(0,*) me,'Realloc blk',isz if (debug) write(0,*) me,'Realloc blk',isz
@ -434,18 +361,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
if (i_ovr < n_ovr) then if (i_ovr < n_ovr) then
if (tot_elem > 1) then if (tot_elem > 1) then
!!$ write(0,*) me,'Realloc temp',tot_elem+2 call imsr(tot_elem,works(idxs+1))
if (tot_elem+2 > size(temp)) then
!!$ write(0,*) me,'Realloc temp',tot_elem+2
deallocate(temp)
allocate(temp(tot_elem+2),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
Call mrgsrt(tot_elem,works(idxs+1),temp,info)
If (info == 0) Call ireordv1(tot_elem,works(idxs+1),temp)
lx = works(idxs+1) lx = works(idxs+1)
i = 1 i = 1
j = 1 j = 1
@ -546,7 +462,18 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
End If End If
desc_p%glob_to_loc(idx)=n_col desc_p%glob_to_loc(idx)=n_col
desc_p%loc_to_glob(n_col)=idx desc_p%loc_to_glob(n_col)=idx
If((counter_t+3) > Size(t_halo_in))Write(0,*)'bingo' If ((counter_t+3) > Size(t_halo_in)) then
isz = max((3*Size(t_halo_in))/2,(counter_t+3+1000))
if (debug) write(0,*) me,'Realloc ovr_el',isz
call psb_realloc(isz,t_halo_in,info)
if (info /= 0) then
info=4010
ch_err='psrealloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end If
t_halo_in(counter_t)=proc_id t_halo_in(counter_t)=proc_id
t_halo_in(counter_t+1)=1 t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col t_halo_in(counter_t+2)=n_col
@ -591,7 +518,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Calling Crea_Halo' if (debug) write(0,*) me,'Calling Crea_Halo'
call psi_crea_index(desc_p,t_halo_in,t_halo_out,.false.,info) call psi_crea_index(desc_p,t_halo_in,t_halo_out,.false.,nxch,nsnd,nrcv,info)
if (debug) then if (debug) then
write(0,*) me,'Done Crea_Index' write(0,*) me,'Done Crea_Index'
@ -611,7 +538,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
tch = tch +(t3-t2) tch = tch +(t3-t2)
End Do End Do
t1 = mpi_wtime() t1 = mpi_wtime()
call FreePairSearchTree(pairtree)
desc_p%matrix_data(psb_m_)=desc_a%matrix_data(psb_m_) desc_p%matrix_data(psb_m_)=desc_a%matrix_data(psb_m_)
desc_p%matrix_data(psb_n_)=desc_a%matrix_data(psb_n_) desc_p%matrix_data(psb_n_)=desc_a%matrix_data(psb_n_)
@ -627,35 +553,12 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
! !
if (debug) then if (debug) then
write(0,*) 'psb_cdasb: converting indexes' write(0,*) 'psb_cdovrbld: converting indexes'
call psb_barrier(ictxt) call psb_barrier(ictxt)
end if end if
!.... convert comunication stuctures.... !.... convert comunication stuctures....
! first the halo index
call psi_crea_index(desc_p,tmp_halo,&
& desc_p%halo_index,.false.,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
! then the overlap index
call psi_crea_index(desc_p,tmp_ovr_idx,&
& desc_p%ovrlap_index,.true.,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_index')
goto 9999
end if
! next is the ovrlap_elem index call psi_cnv_dsc(tmp_halo,tmp_ovr_idx,desc_p,info)
call psi_crea_ovr_elem(desc_p%ovrlap_index,desc_p%ovrlap_elem)
! finally bnd_elem
call psi_crea_bnd_elem(desc_p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psi_crea_bnd_elem')
goto 9999
end if
! Ok, register into MATRIX_DATA & free temporary work areas ! Ok, register into MATRIX_DATA & free temporary work areas
desc_p%matrix_data(psb_dec_type_) = psb_desc_asb_ desc_p%matrix_data(psb_dec_type_) = psb_desc_asb_

@ -54,10 +54,9 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
use psb_tools_mod, only : psb_glob_to_loc, psb_loc_to_glob use psb_tools_mod, only : psb_glob_to_loc, psb_loc_to_glob
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use mpi
Implicit None Implicit None
include 'mpif.h'
Type(psb_zspmat_type),Intent(in) :: a Type(psb_zspmat_type),Intent(in) :: a
Type(psb_zspmat_type),Intent(inout) :: blk Type(psb_zspmat_type),Intent(inout) :: blk
Type(psb_desc_type),Intent(in) :: desc_a Type(psb_desc_type),Intent(in) :: desc_a

Loading…
Cancel
Save