#$=head1 NAME
#$
#$invint2 - Inverse linear interpolation;
#$
#$=head1 SYNOPSIS
#$
#$C<call invint1(niter,coord,ord,o1,d1,mm,mmov,eps,aa,doprec)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item niter - integer     
#$
#$      Number of itterations
#$
#$=item coord -  C<real(:)> 
#$
#$      Coordinates
#$
#$=item o1,d1 -  real       
#$
#$      First model position and sampling
#$
#$=item ord   -  C<real(:)> 
#$
#$      Data values
#$
#$=item mm    -  C<real(:)> 
#$
#$      Output model
#$
#$=item mmov  -  C<real(:,:)> 
#$
#$      Model movie
#$
#$=item eps   -  real         
#$
#$      Epsilon value if doing preconditioning
#$
#$=item aa    -  type(filter)  
#$
#$      Preconditioning operator
#$
#$=item prec  -  logical      
#$
#$      Whether or not to do preconditioning
#$
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Perform inverse linear interpolation 
#$
#$
#$=head1 SEE ALSO
#$
#$L<lint1>,L<helicon>,L<polydiv>,L<solver>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$
module invint2 {   # Inverse linear interpolation
  use lint1
  use helicon                           # regularized    by helix   filtering
  use polydiv                           # preconditioned by inverse filtering
  use cgstep_mod + solver_mod
contains
  subroutine invint1( niter, coord,ord, o1,d1, mm,mmov, eps, aa, doprec) {
    logical,                  intent( in)  :: doprec
    integer,                  intent( in)  :: niter
    real,                     intent( in)  :: o1, d1, eps
    real,    dimension( :),   intent( in)  :: ord  
    type( filter),            intent( in)  :: aa
    real,    dimension( :),   intent( out) :: mm             
    real,    dimension( :,:), intent( out) :: mmov         # model movie
    real,    dimension( :),   pointer      :: coord        # coordinate
    call lint1_init( o1, d1, coord)
    if( doprec) {                                          # preconditioning
       call polydiv_init( size(mm), aa)
       call solver_prec( lint1_lop, cgstep, niter = niter, x = mm, dat = ord,  
                    prec = polydiv_lop, nprec = size(mm), eps = eps, 
                    xmov = mmov)
       call polydiv_close()
    } else {                                               # regularization
       call helicon_init( aa)
       call solver_reg(  lint1_lop, cgstep, niter = niter, x = mm, dat = ord, 
                    reg = helicon_lop, nreg = size(mm), eps = eps,
                    xmov = mmov)
    }                                               
    call cgstep_close()
  }
}