mld2p4-2:

mlprec/impl/Makefile
 mlprec/impl/mld_cumf_interface.c
 mlprec/impl/mld_sumf_interface.c
 mlprec/mld_c_umf_solver.F90
 mlprec/mld_s_umf_solver.F90
 tests/pdegen/runs/mld_pde2d.inp

Get rid of nonsensical single-precision interfaces to UMFPACK.
stopcriterion
Salvatore Filippone 8 years ago
parent aa5703e62a
commit c99c587c44

@ -61,9 +61,9 @@ OUTEROBJS=$(SOUTEROBJS) $(DOUTEROBJS) $(COUTEROBJS) $(ZOUTEROBJS)
F90OBJS=$(OUTEROBJS) $(INNEROBJS)
COBJS= mld_sslu_interface.o mld_sumf_interface.o \
COBJS= mld_sslu_interface.o \
mld_dslu_interface.o mld_dumf_interface.o \
mld_cslu_interface.o mld_cumf_interface.o \
mld_cslu_interface.o \
mld_zslu_interface.o mld_zumf_interface.o
OBJS=$(F90OBJS) $(COBJS) $(MPCOBJS)

@ -1,195 +0,0 @@
/*
*
* MLD2P4 version 2.0
* MultiLevel Domain Decomposition Parallel Preconditioners Package
* based on PSBLAS (Parallel Sparse BLAS version 3.3)
*
* (C) Copyright 2008, 2010, 2012, 2015
*
* Salvatore Filippone University of Rome Tor Vergata
* Alfredo Buttari CNRS-IRIT, Toulouse
* Pasqua D'Ambra ICAR-CNR, Naples
* Daniela di Serafino Second University of Naples
*
* 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 MLD2P4 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 MLD2P4 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.
*
*
* File: mld_cumf_interface.c
*
* Functions: mld_cumf_fact_, mld_cumf_solve_, mld_cumf_free_.
*
* This file is an interface to the UMFPACK routines for sparse factorization and
* solve. It was obtained by adapting umfpack_zi_demo under the original UMFPACK
* copyright terms reproduced below.
*
*/
/* =====================
UMFPACK Version 4.4 (Jan. 28, 2005), Copyright (c) 2005 by Timothy A.
Davis. All Rights Reserved.
UMFPACK License:
Your use or distribution of UMFPACK or any modified version of
UMFPACK implies that you agree to this License.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use or copy this program, provided
that the Copyright, this License, and the Availability of the original
version is retained on all copies. User documentation of any code that
uses UMFPACK or any modified version of UMFPACK code must cite the
Copyright, this License, the Availability note, and "Used by permission."
Permission to modify the code and to distribute modified code is granted,
provided the Copyright, this License, and the Availability note are
retained, and a notice that the code was modified is included. This
software was developed with support from the National Science Foundation,
and is provided to you free of charge.
Availability:
http://www.cise.ufl.edu/research/sparse/umfpack
*/
#include <stdio.h>
#ifdef Have_UMF_
#include "umfpack.h"
#endif
int mld_cumf_fact(int n, int nnz,
float *values, int *rowind, int *colptr,
void **symptr, void **numptr,
long long int *ssize,
long long int *nsize)
{
#ifdef Have_UMF_
double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL];
void *Symbolic, *Numeric ;
int i, info;
umfpack_si_defaults(Control);
info = umfpack_ci_symbolic (n, n, colptr, rowind, values, NULL, &Symbolic,
Control, Info);
if ( info == UMFPACK_OK ) {
info = 0;
} else {
printf("umfpack_ci_symbolic() error returns INFO= %d\n", info);
*symptr = (void *) NULL;
*numptr = (void *) NULL;
return -11;
}
*symptr = Symbolic;
*ssize = Info[UMFPACK_SYMBOLIC_SIZE];
*ssize *= Info[UMFPACK_SIZE_OF_UNIT];
info = umfpack_ci_numeric (colptr, rowind, values, NULL, Symbolic, &Numeric,
Control, Info) ;
if ( info == UMFPACK_OK ) {
info = 0;
*numptr = Numeric;
*nsize = Info[UMFPACK_NUMERIC_SIZE];
*nsize *= Info[UMFPACK_SIZE_OF_UNIT];
} else {
printf("umfpack_ci_numeric() error returns INFO= %d\n", info);
info = -12;
*numptr = NULL;
}
return info;
#else
fprintf(stderr," UMF Not available for single precision.\n");
return -1;
#endif
}
int mld_cumf_solve(int itrans, int n,
float *x, float *b, int ldb,
void *numptr)
{
#ifdef Have_UMF_
double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL];
void *Symbolic, *Numeric ;
int i,trans, info;
umfpack_si_defaults(Control);
Control[UMFPACK_IRSTEP]=0;
if (itrans == 0) {
trans = UMFPACK_A;
} else if (itrans ==1) {
trans = UMFPACK_At;
} else {
trans = UMFPACK_A;
}
info = umfpack_ci_solve(trans,NULL,NULL,NULL,NULL,
x,NULL,b,NULL, numptr,Control,Info);
return info;
#else
fprintf(stderr," UMF Not available for single precision.\n");
return -1;
#endif
}
int mld_cumf_free(void *symptr, void *numptr)
{
#ifdef Have_UMF_
void *Symbolic, *Numeric ;
Symbolic = symptr;
Numeric = numptr;
umfpack_ci_free_numeric(&Numeric);
umfpack_ci_free_symbolic(&Symbolic);
return 0;
#else
fprintf(stderr," UMF Not available for single precision.\n");
return -1;
#endif
}

@ -1,201 +0,0 @@
/*
*
* MLD2P4 version 2.0
* MultiLevel Domain Decomposition Parallel Preconditioners Package
* based on PSBLAS (Parallel Sparse BLAS version 3.3)
*
* (C) Copyright 2008, 2010, 2012, 2015
*
* Salvatore Filippone University of Rome Tor Vergata
* Alfredo Buttari CNRS-IRIT, Toulouse
* Pasqua D'Ambra ICAR-CNR, Naples
* Daniela di Serafino Second University of Naples
*
* 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 MLD2P4 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 MLD2P4 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.
*
*
* File: mld_umf_interface.c
*
* Functions: mld_sumf_fact_, mld_sumf_solve_, mld_umf_free_.
*
* This file is an interface to the UMFPACK routines for sparse factorization and
* solve. It was obtained by adapting umfpack_di_demo under the original UMFPACK
* copyright terms reproduced below.
*
*/
/* =====================
UMFPACK Version 4.4 (Jan. 28, 2005), Copyright (c) 2005 by Timothy A.
Davis. All Rights Reserved.
UMFPACK License:
Your use or distribution of UMFPACK or any modified version of
UMFPACK implies that you agree to this License.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use or copy this program, provided
that the Copyright, this License, and the Availability of the original
version is retained on all copies. User documentation of any code that
uses UMFPACK or any modified version of UMFPACK code must cite the
Copyright, this License, the Availability note, and "Used by permission."
Permission to modify the code and to distribute modified code is granted,
provided the Copyright, this License, and the Availability note are
retained, and a notice that the code was modified is included. This
software was developed with support from the National Science Foundation,
and is provided to you free of charge.
Availability:
http://www.cise.ufl.edu/research/sparse/umfpack
*/
#include <stdio.h>
/* Currently no single precision version in UMFPACK */
#ifdef Have_UMF_
#undef Have_UMF_
#endif
#ifdef Have_UMF_
#include "umfpack.h"
#endif
int mld_sumf_fact(int n, int nnz,
float *values, int *rowind, int *colptr,
void *symptr,
void *numptr,
long long int *ssize,
long long int *nsize)
{
#ifdef Have_UMF_
float Info [UMFPACK_INFO], Control [UMFPACK_CONTROL];
void *Symbolic, *Numeric ;
int i,info ;
umfpack_si_defaults(Control);
info = umfpack_si_symbolic (n, n, colptr, rowind, values, &Symbolic,
Control, Info);
if ( info == UMFPACK_OK ) {
info = 0;
} else {
printf("umfpack_si_symbolic() error returns INFO= %d\n", info);
*symptr = (void *) NULL;
*numptr = (void *) NULL;
return -11;
}
*symptr = Symbolic;
*ssize = Info[UMFPACK_SYMBOLIC_SIZE];
*ssize *= Info[UMFPACK_SIZE_OF_UNIT];
info = umfpack_si_numeric (colptr, rowind, values, Symbolic, &Numeric,
Control, Info) ;
if ( info == UMFPACK_OK ) {
info = 0;
*numptr = Numeric;
*nsize = Info[UMFPACK_NUMERIC_SIZE];
*nsize *= Info[UMFPACK_SIZE_OF_UNIT];
} else {
printf("umfpack_si_numeric() error returns INFO= %d\n", info);
info = -12;
*numptr = NULL;
}
return info;
#else
fprintf(stderr," UMF Not available for single precision.\n");
return -1;
#endif
}
int mld_sumf_solve(int itrans, int n,
float *x, float *b, int ldb,
void *numptr)
{
#ifdef Have_UMF_
float Info [UMFPACK_INFO], Control [UMFPACK_CONTROL];
void *Symbolic, *Numeric ;
int i,trans, info;
umfpack_si_defaults(Control);
Control[UMFPACK_IRSTEP]=0;
if (itrans == 0) {
trans = UMFPACK_A;
} else if (itrans ==1) {
trans = UMFPACK_At;
} else {
trans = UMFPACK_A;
}
info = umfpack_si_solve(trans,NULL,NULL,NULL,
x,b,(void *) *numptr,Control,Info);
return info;
#else
fprintf(stderr," UMF Not available for single precision.\n");
return -1;
#endif
}
int mld_sumf_free(void *symptr, void *numptr)
{
#ifdef Have_UMF_
void *Symbolic, *Numeric ;
Symbolic = symptr;
Numeric = numptr;
if (numptr != NULL) umfpack_si_free_numeric(&Numeric);
if (symptr != NULL) umfpack_si_free_symbolic(&Symbolic);
return 0;
#else
fprintf(stderr," UMF Not available for single precision.\n");
return -1;
#endif
}

@ -1,426 +0,0 @@
!!$
!!$
!!$ MLD2P4 version 2.1
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.4)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015, 2017
!!$
!!$ Salvatore Filippone Cranfield University
!!$ Ambra Abdullahi Hassan University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ 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 MLD2P4 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 MLD2P4 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.
!!$
!!$
!
! File: mld_c_umf_solver_mod.f90
!
! Module: mld_c_umf_solver_mod
!
! This module defines:
! - the mld_c_umf_solver_type data structure containing the ingredients
! to interface with the UMFPACK package.
! 1. The factorization is restricted to the diagonal block of the
! current image.
!
module mld_c_umf_solver
use iso_c_binding
use mld_c_base_solver_mod
#if defined(LONG_INTEGERS)
type, extends(mld_c_base_solver_type) :: mld_c_umf_solver_type
end type mld_c_umf_solver_type
#else
type, extends(mld_c_base_solver_type) :: mld_c_umf_solver_type
type(c_ptr) :: symbolic=c_null_ptr, numeric=c_null_ptr
integer(c_long_long) :: symbsize=0, numsize=0
contains
procedure, pass(sv) :: build => c_umf_solver_bld
procedure, pass(sv) :: apply_a => c_umf_solver_apply
procedure, pass(sv) :: apply_v => c_umf_solver_apply_vect
procedure, pass(sv) :: free => c_umf_solver_free
procedure, pass(sv) :: descr => c_umf_solver_descr
procedure, pass(sv) :: sizeof => c_umf_solver_sizeof
#if defined(HAVE_FINAL)
final :: c_umf_solver_finalize
#endif
end type mld_c_umf_solver_type
private :: c_umf_solver_bld, c_umf_solver_apply, &
& c_umf_solver_free, c_umf_solver_descr, &
& c_umf_solver_sizeof, c_umf_solver_apply_vect
#if defined(HAVE_FINAL)
private :: c_umf_solver_finalize
#endif
interface
function mld_cumf_fact(n,nnz,values,rowind,colptr,&
& symptr,numptr,ssize,nsize)&
& bind(c,name='mld_cumf_fact') result(info)
use iso_c_binding
integer(c_int), value :: n,nnz
integer(c_int) :: info
integer(c_long_long) :: ssize, nsize
integer(c_int) :: rowind(*),colptr(*)
complex(c_float_complex) :: values(*)
type(c_ptr) :: symptr, numptr
end function mld_cumf_fact
end interface
interface
function mld_cumf_solve(itrans,n,x, b, ldb, numptr)&
& bind(c,name='mld_cumf_solve') result(info)
use iso_c_binding
integer(c_int) :: info
integer(c_int), value :: itrans,n,ldb
complex(c_float_complex) :: x(*), b(ldb,*)
type(c_ptr), value :: numptr
end function mld_cumf_solve
end interface
interface
function mld_cumf_free(symptr, numptr)&
& bind(c,name='mld_cumf_free') result(info)
use iso_c_binding
integer(c_int) :: info
type(c_ptr), value :: symptr, numptr
end function mld_cumf_free
end interface
contains
subroutine c_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
& trans,work,info,init,initu)
use psb_base_mod
implicit none
type(psb_desc_type), intent(in) :: desc_data
class(mld_c_umf_solver_type), intent(inout) :: sv
complex(psb_spk_),intent(inout) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
complex(psb_spk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
character, intent(in), optional :: init
complex(psb_spk_),intent(inout), optional :: initu(:)
integer :: n_row,n_col
complex(psb_spk_), pointer :: ww(:)
integer :: ictxt,np,me,i, err_act
character :: trans_
character(len=20) :: name='c_umf_solver_apply'
call psb_erractionsave(err_act)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
!
! For non-iterative solvers, init and initu are ignored.
!
n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then
ww => work(1:n_col)
else
allocate(ww(n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
endif
select case(trans_)
case('N')
info = mld_cumf_solve(0,n_row,ww,x,n_row,sv%numeric)
case('T')
!
! Note: with UMF, 1 meand Ctranspose, 2 means transpose
! even for complex data.
!
if (psb_c_is_complex_) then
info = mld_cumf_solve(2,n_row,ww,x,n_row,sv%numeric)
else
info = mld_cumf_solve(1,n_row,ww,x,n_row,sv%numeric)
end if
case('C')
info = mld_cumf_solve(1,n_row,ww,x,n_row,sv%numeric)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
goto 9999
endif
if (n_col > size(work)) then
deallocate(ww)
endif
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine c_umf_solver_apply
subroutine c_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
& trans,work,info,init,initu)
use psb_base_mod
implicit none
type(psb_desc_type), intent(in) :: desc_data
class(mld_c_umf_solver_type), intent(inout) :: sv
type(psb_c_vect_type),intent(inout) :: x
type(psb_c_vect_type),intent(inout) :: y
complex(psb_spk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu
integer :: err_act
character(len=20) :: name='c_umf_solver_apply_vect'
call psb_erractionsave(err_act)
info = psb_success_
!
! For non-iterative solvers, init and initu are ignored.
!
call x%v%sync()
call y%v%sync()
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
call y%v%set_host()
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine c_umf_solver_apply_vect
subroutine c_umf_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold)
use psb_base_mod
Implicit None
! Arguments
type(psb_cspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(mld_c_umf_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
type(psb_cspmat_type), intent(in), target, optional :: b
class(psb_c_base_sparse_mat), intent(in), optional :: amold
class(psb_c_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
type(psb_cspmat_type) :: atmp
type(psb_c_csc_sparse_mat) :: acsc
integer :: n_row,n_col, nrow_a, nztota
integer :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='c_umf_solver_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
if (psb_toupper(upd) == 'F') then
call a%cscnv(atmp,info,type='coo')
call psb_rwextd(n_row,atmp,info,b=b)
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
call atmp%mv_to(acsc)
nrow_a = acsc%get_nrows()
nztota = acsc%get_nzeros()
! Fix the entres to call C-base UMFPACK.
acsc%ia(:) = acsc%ia(:) - 1
acsc%icp(:) = acsc%icp(:) - 1
info = mld_cumf_fact(nrow_a,nztota,acsc%val,&
& acsc%ia,acsc%icp,sv%symbolic,sv%numeric,&
& sv%symbsize,sv%numsize)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='mld_cumf_fact'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call acsc%free()
call atmp%free()
else
! ?
info=psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine c_umf_solver_bld
subroutine c_umf_solver_free(sv,info)
Implicit None
! Arguments
class(mld_c_umf_solver_type), intent(inout) :: sv
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='c_umf_solver_free'
call psb_erractionsave(err_act)
info = mld_cumf_free(sv%symbolic,sv%numeric)
if (info /= psb_success_) goto 9999
sv%symbolic = c_null_ptr
sv%numeric = c_null_ptr
sv%symbsize = 0
sv%numsize = 0
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine c_umf_solver_free
#if defined(HAVE_FINAL)
subroutine c_umf_solver_finalize(sv)
Implicit None
! Arguments
type(mld_c_umf_solver_type), intent(inout) :: sv
integer :: info
Integer :: err_act
character(len=20) :: name='c_umf_solver_finalize'
call sv%free(info)
return
end subroutine c_umf_solver_finalize
#endif
subroutine c_umf_solver_descr(sv,info,iout,coarse)
Implicit None
! Arguments
class(mld_c_umf_solver_type), intent(in) :: sv
integer, intent(out) :: info
integer, intent(in), optional :: iout
logical, intent(in), optional :: coarse
! Local variables
integer :: err_act
integer :: ictxt, me, np
character(len=20), parameter :: name='mld_c_umf_solver_descr'
integer :: iout_
call psb_erractionsave(err_act)
info = psb_success_
if (present(iout)) then
iout_ = iout
else
iout_ = 6
endif
write(iout_,*) ' UMFPACK Sparse Factorization Solver. '
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine c_umf_solver_descr
function c_umf_solver_sizeof(sv) result(val)
implicit none
! Arguments
class(mld_c_umf_solver_type), intent(in) :: sv
integer(psb_long_int_k_) :: val
integer :: i
val = 2*psb_sizeof_long_int
val = val + sv%symbsize
val = val + sv%numsize
return
end function c_umf_solver_sizeof
#endif
end module mld_c_umf_solver

@ -1,426 +0,0 @@
!!$
!!$
!!$ MLD2P4 version 2.1
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.4)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015, 2017
!!$
!!$ Salvatore Filippone Cranfield University
!!$ Ambra Abdullahi Hassan University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ 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 MLD2P4 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 MLD2P4 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.
!!$
!!$
!
! File: mld_s_umf_solver_mod.f90
!
! Module: mld_s_umf_solver_mod
!
! This module defines:
! - the mld_s_umf_solver_type data structure containing the ingredients
! to interface with the UMFPACK package.
! 1. The factorization is restricted to the diagonal block of the
! current image.
!
module mld_s_umf_solver
use iso_c_binding
use mld_s_base_solver_mod
#if defined(LONG_INTEGERS)
type, extends(mld_s_base_solver_type) :: mld_s_umf_solver_type
end type mld_s_umf_solver_type
#else
type, extends(mld_s_base_solver_type) :: mld_s_umf_solver_type
type(c_ptr) :: symbolic=c_null_ptr, numeric=c_null_ptr
integer(c_long_long) :: symbsize=0, numsize=0
contains
procedure, pass(sv) :: build => s_umf_solver_bld
procedure, pass(sv) :: apply_a => s_umf_solver_apply
procedure, pass(sv) :: apply_v => s_umf_solver_apply_vect
procedure, pass(sv) :: free => s_umf_solver_free
procedure, pass(sv) :: descr => s_umf_solver_descr
procedure, pass(sv) :: sizeof => s_umf_solver_sizeof
#if defined(HAVE_FINAL)
final :: s_umf_solver_finalize
#endif
end type mld_s_umf_solver_type
private :: s_umf_solver_bld, s_umf_solver_apply, &
& s_umf_solver_free, s_umf_solver_descr, &
& s_umf_solver_sizeof, s_umf_solver_apply_vect
#if defined(HAVE_FINAL)
private :: s_umf_solver_finalize
#endif
interface
function mld_sumf_fact(n,nnz,values,rowind,colptr,&
& symptr,numptr,ssize,nsize)&
& bind(c,name='mld_sumf_fact') result(info)
use iso_c_binding
integer(c_int), value :: n,nnz
integer(c_int) :: info
integer(c_long_long) :: ssize, nsize
integer(c_int) :: rowind(*),colptr(*)
real(c_float) :: values(*)
type(c_ptr) :: symptr, numptr
end function mld_sumf_fact
end interface
interface
function mld_sumf_solve(itrans,n,x, b, ldb, numptr)&
& bind(c,name='mld_sumf_solve') result(info)
use iso_c_binding
integer(c_int) :: info
integer(c_int), value :: itrans,n,ldb
real(c_float) :: x(*), b(ldb,*)
type(c_ptr), value :: numptr
end function mld_sumf_solve
end interface
interface
function mld_sumf_free(symptr, numptr)&
& bind(c,name='mld_sumf_free') result(info)
use iso_c_binding
integer(c_int) :: info
type(c_ptr), value :: symptr, numptr
end function mld_sumf_free
end interface
contains
subroutine s_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
& trans,work,info,init,initu)
use psb_base_mod
implicit none
type(psb_desc_type), intent(in) :: desc_data
class(mld_s_umf_solver_type), intent(inout) :: sv
real(psb_spk_),intent(inout) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
real(psb_spk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
character, intent(in), optional :: init
real(psb_spk_),intent(inout), optional :: initu(:)
integer :: n_row,n_col
real(psb_spk_), pointer :: ww(:)
integer :: ictxt,np,me,i, err_act
character :: trans_
character(len=20) :: name='s_umf_solver_apply'
call psb_erractionsave(err_act)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
!
! For non-iterative solvers, init and initu are ignored.
!
n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then
ww => work(1:n_col)
else
allocate(ww(n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
& a_err='real(psb_spk_)')
goto 9999
end if
endif
select case(trans_)
case('N')
info = mld_sumf_solve(0,n_row,ww,x,n_row,sv%numeric)
case('T')
!
! Note: with UMF, 1 meand Ctranspose, 2 means transpose
! even for complex data.
!
if (psb_s_is_complex_) then
info = mld_sumf_solve(2,n_row,ww,x,n_row,sv%numeric)
else
info = mld_sumf_solve(1,n_row,ww,x,n_row,sv%numeric)
end if
case('C')
info = mld_sumf_solve(1,n_row,ww,x,n_row,sv%numeric)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
goto 9999
endif
if (n_col > size(work)) then
deallocate(ww)
endif
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine s_umf_solver_apply
subroutine s_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
& trans,work,info,init,initu)
use psb_base_mod
implicit none
type(psb_desc_type), intent(in) :: desc_data
class(mld_s_umf_solver_type), intent(inout) :: sv
type(psb_s_vect_type),intent(inout) :: x
type(psb_s_vect_type),intent(inout) :: y
real(psb_spk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu
integer :: err_act
character(len=20) :: name='s_umf_solver_apply_vect'
call psb_erractionsave(err_act)
info = psb_success_
!
! For non-iterative solvers, init and initu are ignored.
!
call x%v%sync()
call y%v%sync()
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
call y%v%set_host()
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine s_umf_solver_apply_vect
subroutine s_umf_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold)
use psb_base_mod
Implicit None
! Arguments
type(psb_sspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(mld_s_umf_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
type(psb_sspmat_type), intent(in), target, optional :: b
class(psb_s_base_sparse_mat), intent(in), optional :: amold
class(psb_s_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
type(psb_sspmat_type) :: atmp
type(psb_s_csc_sparse_mat) :: acsc
integer :: n_row,n_col, nrow_a, nztota
integer :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='s_umf_solver_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
if (psb_toupper(upd) == 'F') then
call a%cscnv(atmp,info,type='coo')
call psb_rwextd(n_row,atmp,info,b=b)
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
call atmp%mv_to(acsc)
nrow_a = acsc%get_nrows()
nztota = acsc%get_nzeros()
! Fix the entres to call C-base UMFPACK.
acsc%ia(:) = acsc%ia(:) - 1
acsc%icp(:) = acsc%icp(:) - 1
info = mld_sumf_fact(nrow_a,nztota,acsc%val,&
& acsc%ia,acsc%icp,sv%symbolic,sv%numeric,&
& sv%symbsize,sv%numsize)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='mld_sumf_fact'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call acsc%free()
call atmp%free()
else
! ?
info=psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine s_umf_solver_bld
subroutine s_umf_solver_free(sv,info)
Implicit None
! Arguments
class(mld_s_umf_solver_type), intent(inout) :: sv
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='s_umf_solver_free'
call psb_erractionsave(err_act)
info = mld_sumf_free(sv%symbolic,sv%numeric)
if (info /= psb_success_) goto 9999
sv%symbolic = c_null_ptr
sv%numeric = c_null_ptr
sv%symbsize = 0
sv%numsize = 0
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine s_umf_solver_free
#if defined(HAVE_FINAL)
subroutine s_umf_solver_finalize(sv)
Implicit None
! Arguments
type(mld_s_umf_solver_type), intent(inout) :: sv
integer :: info
Integer :: err_act
character(len=20) :: name='s_umf_solver_finalize'
call sv%free(info)
return
end subroutine s_umf_solver_finalize
#endif
subroutine s_umf_solver_descr(sv,info,iout,coarse)
Implicit None
! Arguments
class(mld_s_umf_solver_type), intent(in) :: sv
integer, intent(out) :: info
integer, intent(in), optional :: iout
logical, intent(in), optional :: coarse
! Local variables
integer :: err_act
integer :: ictxt, me, np
character(len=20), parameter :: name='mld_s_umf_solver_descr'
integer :: iout_
call psb_erractionsave(err_act)
info = psb_success_
if (present(iout)) then
iout_ = iout
else
iout_ = 6
endif
write(iout_,*) ' UMFPACK Sparse Factorization Solver. '
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine s_umf_solver_descr
function s_umf_solver_sizeof(sv) result(val)
implicit none
! Arguments
class(mld_s_umf_solver_type), intent(in) :: sv
integer(psb_long_int_k_) :: val
integer :: i
val = 2*psb_sizeof_long_int
val = val + sv%symbsize
val = val + sv%numsize
return
end function s_umf_solver_sizeof
#endif
end module mld_s_umf_solver

@ -13,7 +13,7 @@ ML ! Preconditioner NONE JACOBI BJAC AS ML
-4 ! If ML: Prescribed number of levels; if <0, ignore it and use coarse size.
-8000 ! If ML: Target coarse size. If <0, then use library default
-1.5d0 ! If ML: Minimum aggregation ratio; if <0 use library default
-0.10d0 ! If ML: Smoother Aggregation Threshold: >= 0.0 default if <0
-0.08d0 ! If ML: Smoother Aggregation Threshold: >= 0.0 default if <0
-20 ! If ML: Maximum acceptable number of levels; if <0 use library default
SMOOTHED ! Type of aggregation: SMOOTHED, UNSMOOTHED, MINENERGY
DEC ! Type of aggregation: DEC SYMDEC

Loading…
Cancel
Save