text:
c
c     ****************************************************************
c     *                                                              *
c     *  assemble & solve linear equations for a Newton iteration    *
c     *                                                              *
c     *                       written by  : rhd                      *
c     *                   last modified : 11/19/2019 rhd             *
c     *                                                              *
c     ****************************************************************
c
      subroutine drive_assemble_solve( first_solve, now_iteration,
     &                                 suggested_new_precond )
      use global_data ! old common.main
c
      use elem_block_data, only : edest_blocks
      use main_data, only : repeat_incid, modified_mpcs,
     &                      asymmetric_assembly, force_solver_rebuild
      use stiffness_data, only : ncoeff, big_ncoeff, k_coeffs,
     &                           k_indexes,
     &                           ncoeff_from_assembled_profile
      use mod_mpc, only : tied_con_mpcs_constructed, mpcs_exist
      use hypre_parameters, only: precond_fail_count, hyp_trigger_step
      use performance_data
      use distributed_stiffness_data, only: parallel_assembly_used,
     &                                   distributed_stiffness_used
c
      implicit none
c
c                    parameter declarations
c
      integer :: now_iteration
      logical :: first_solve, suggested_new_precond
c
c                    locals
c
      integer :: i, k, error_code, nnz, mkl_threads, num_enode_dof,
     &           num_struct_dof, iresult, itype, num_terms, i_offset,
     &    eqn_num
      integer, save :: neqns, old_neqns, old_ncoeff
      integer, external :: curr_neqns
      integer :: code_vec(mxedof,mxvl), edest(mxedof,mxconn)
      integer, allocatable :: k_ptrs(:)
      integer, allocatable, save :: dof_eqn_map(:), eqn_node_map(:),
     &                              save_k_indexes(:), save_k_ptrs(:)

      logical :: new_size
      logical, save :: cpu_stats, save_solver
      logical, parameter :: local_debug = .false.,
     &     local_debug2 = .false., local_debug3 = .false.
c
      real, external :: wcputime
      double precision, parameter :: zero = 0.0d00
      double precision, allocatable :: p_vec(:), k_diag(:), u_vec(:)
c
      data old_neqns, old_ncoeff, cpu_stats, save_solver
     &     / 0, 0, .true., .false. /
c
      if( local_debug ) write(*,*) '... drive_assem_solve ... @ 1'
      if( .not. show_details ) cpu_stats = .false.
c
      num_enode_dof  = iprops(4,1) ! always = 3 in warp3d
      num_struct_dof = nonode * num_enode_dof
      new_size       = .false.
c
c              new capability for any part of code to force
c              reconstruction of all data structures for the solver.
c              E.g. releasing of MPCs where this chance cannot
c              be detected -- essentially a massive data reconstruction.
c
      if( force_solver_rebuild ) then
          new_size = .true.
          force_solver_rebuild = .false.
          cpu_stats = .true.
      end if
c
c              stiffness data structure development. can skip
c              if not 1st solution for this load(time) step.
c              can also skip if # eqns and sparsity unchanged
c              from prior load(time) step.

      if( first_solve ) then ! 1st for this load(time) step
c
c           1. build sparsity data structures
c              ------------------------------
c
c              based on decisions/checks made in mnralg,
c              setup for threaded or MPI parallel assembly by
c              generating the sparsity data structures
c
        call t_start_assembly( start_assembly_step ) ! timing
        if( .not. parallel_assembly_used ) then
          call ds_setup_sparsity_local( iresult, use_mpi, numprocs )
          if( iresult .eq. 1 ) then ! no equations to solve
             idu(1:nodof) = zero
             return
          end if
        end if
        if( local_debug ) write(*,*) '... drive_assem_solve  @ 2'
        if( parallel_assembly_used ) then
          call ds_setup_sparsity_distributed( iresult )
          if( iresult .eq. 0 ) then  ! no equations to solve
             idu(1:nodof) = zero
             return
          end if
        end if
        call t_end_assembly( assembly_total, start_assembly_step )
        if( local_debug ) write(*,*) '... drive_assem_solve  @ 3'

c
      end if   ! on first_solve
c
c          2.  assemble the equations and rh side here
c              ---------------------------------------
c
c             => threads only assembly
c
c                - symmetric assembly with sparse (MKL)
c                  direct/iterative solvers. Can have MPC
c                  from user or tied-contact.
c
c                - asymmetric assembly with sparse (MKL)
c                  direct/iterative solvers. The sparsity is symmetric
c                  but not coefficients, e.g., from crystal plasticity
c                  material model
c
c                - a symmetric/asymmetric assembly on root followed by
c                  distribution to MPI ranks for solution by hypre
c
c                or
c
c                - a symmetric/asymmetric assembly on root follwed
c                  by use of MPI cpardiso direct solver on ranks
c
c             => distributed assembly across MPI ranks followed by
c                hypre or CPardiso solve.
c
      call t_start_assembly( start_assembly_step )
      ntimes_assembly = ntimes_assembly + 1
      if( .not. parallel_assembly_used ) then
        call ds_root_assembly
        itype = 2
        if( new_size ) itype = 1
       end if
      if( parallel_assembly_used ) call ds_distributed_assembly
c
c              distributed assembly across MPI ranks followed by
c              do we need a distributed stiffness -- but did not use
c              parallel assembly?  If so, take care of that now.
c              The "true" flag indicates full (not LT) assembly and the
c             'blockrow' parameter means agglomerate by rows and map
c              by blocks of equal nnz size.
c
      if( .not. parallel_assembly_used
     &    .and. distributed_stiffness_used ) then
            call wmpi_alert_slaves(39)
            call distribute_from_assembled(neqns, ncoeff,
     &                  k_diag, p_vec, u_vec, k_coeffs, k_ptrs,
     &                  k_indexes, .true., "blockrow", itype, out )
      end if
c
      if( local_debug ) write(*,*) '... drive_assem_solve @ 4'
      call t_end_assembly( assembly_total, start_assembly_step )
c
c
c          3.  solve the equations
c              -------------------
c
c              computed u_vec has the corrections to du.
c              du is current estimate of total displacement
c              change from n -> n+1
c
c              p_vec was set above and is the set of residual
c              nodal forces at n+1. they include MPC effects as
c              needed for symmetric MKL solvers.
c
c               hypre_solver        => solver_flag .eq. 9
c               pardiso_asymmetric  => solver_flag .eq. 8
c               pardiso_symmetric   => solver_flag .eq. 7
c               cpardiso_symmetric  => solver_flag .eq. 10
c               cpardiso_asymmetric => solver_flag .eq. 11
c
      select case( solver_flag )

      case( 7 ) ! symmetric pardiso
c
        if( asymmetric_assembly ) then
          write(out,9120); call die_gracefully
        end if
        if( local_debug ) write(*,*) '... drive_assem_solve  @ 5a'
        call pardiso_symmetric( neqns, ncoeff, k_diag, p_vec,
     &                          u_vec, k_coeffs, k_ptrs, k_indexes,
     &                          cpu_stats, itype, out,
     &                          solver_out_of_core, solver_memory,
     &                          solver_scr_dir, solver_mkl_iterative )
        if( local_debug ) write(*,*) '... drive_assem_solve @ 5b'

      case( 8 ) ! asymmetric pardiso
c
        if( .not. asymmetric_assembly ) then
          write(out,9110); call die_gracefully
        end if
        if( local_debug ) write(*,*) '... drive_assem_solve  @ 5c'
        call pardiso_unsymmetric( neqns, nnz, k_ptrs, k_indexes,
     &            k_coeffs,  p_vec, u_vec, cpu_stats, itype, out,
     &            solver_mkl_iterative )
        if( local_debug ) write(*,*) '... drive_assem_solve @ 5d'
c

      case( 9 ) ! hypre_solver
c
        call ds_drive_hypre
        if( hyp_trigger_step ) return ! failed and wants Newton
c
      case( 10 ) ! symmetric cpardiso
c
        if( asymmetric_assembly ) then
          write(out,9120); call die_gracefully
        end if
        if( local_debug ) write(*,*) '... drive_assem_solve @ 5e'
        call wmpi_alert_slaves( 3 )
        call cpardiso_symmetric( neqns, ncoeff, k_diag, p_vec,
     &                          u_vec, k_coeffs, k_ptrs, k_indexes,
     &                          cpu_stats, itype, out, myid )
        if( local_debug ) write(*,*) '... drive_assem_solve @ 5f'
c
      case( 11 ) ! asymmetric cpardiso
c
        if( .not. asymmetric_assembly ) then
             write(out,9125); call die_gracefully
        end if
        if( local_debug ) write(*,*) '... drive_assem_solve @ 5g'
        call wmpi_alert_slaves( 31 )
        call cpardiso_unsymmetric( neqns, nnz, k_ptrs, k_indexes,
     &            k_coeffs,  p_vec, u_vec, cpu_stats, itype, out,
     &            myid )
        if( local_debug ) write(*,*) '... drive_assem_solve @ 5h'
c
      case default ! bad solver type
c
            write(out,9301) solver_flag
            call die_gracefully
c
      end select
c
      if( local_debug ) write(out,*) ' @ 6'
      if( .not. parallel_assembly_used )
     &  deallocate( k_diag, k_coeffs, k_ptrs, k_indexes, p_vec )
c
c          4.  for distributed assembly/solve, reorder solution vector
c              -------------------------------------------------------
c
c              to eliminate the effects of the remapping
c              to improve load balancing
c
      call t_start_assembly( start_assembly_step )
      call reorder_soln_vec( u_vec )
      call t_end_assembly( assembly_total, start_assembly_step )
c
c          5.  apply mpc equations
c              -------------------
c
c              if they exist to compute the corrective displacements
c              to du at dependent dofs.
c
c              compute corrections estimate of increment of
c              lagrange multipliers over the step
c
c              these are nodal forces at all MPC dof to enforce the
c              multi-points constraints
c
c              only for symmetric systems. put into their
c              positions in u_vec.
c
      if( tied_con_mpcs_constructed .or. mpcs_exist ) then
         if( local_debug ) write(out,*) ' @ 7'
         call mpcs_apply( u_vec, neqns, nodof, cstmap )
      end if
c
c
c          6.  expand compressed u_vec to structure nodof
c              ------------------------------------------
c
c              u_vec from solver is the reduced number of equations
c              from absolute constraints and MPCs. expand to
c              nodof size for structure.
c
c              Again - these are the corrective displacements to
c              be added to du.
c
      do i = 1, nodof
       if( dof_eqn_map(i) .eq. 0 ) then
          idu(i) = zero
       else
          idu(i) = u_vec(dof_eqn_map(i))
       end if
      end do
c
      if( local_debug2 ) write(out,9900) (i, u_vec(i), i=1, neqns)
c
      deallocate( u_vec ) ! scratch vec used by solvers
c
c          7.  re-define number of non-zero coefficients in the
c              equations as solved.
c              -------------------------------------------------
c
c              needed only if we inserted mpcs. The next time thru
c              the sover this now "correct" count of the number of
c              terms. saves a lot of work during re-insertion of the
c              mpc coefficients.
c
      if( tied_con_mpcs_constructed .or. mpcs_exist )
     &    ncoeff = big_ncoeff
c
      cpu_stats = .false.   ! only print them 1st time thru solver
c
      return
c
 9100 format('>> constrained load vector: ',
     & 1000(/3x,8f12.6) )
 9110 format(1x,'>> FATAL ERROR: Job Aborted.',
     & /,5x,'asymmetric solver requested for symmetric equations')
 9115 format(1x,'>> FATAL ERROR: Job Aborted.',
     & /,5x,'asymmetric solver requested for symmetric equations')
 9120 format(1x,'>> FATAL ERROR: Job Aborted.',
     & /,5x,'symmetric solver requested for asymmetric equations')
 9125 format(1x,'>> FATAL ERROR: Job Aborted.',
     & /,5x,'asymmetric solver requested w/o asymmetric assembly')
 9301 format(1x,
     &'>> FATAL ERROR: Job Aborted.',
     &5x,' Requested equation solver cannot be used',
     &/5x,'Solver_flag is set to:',i12 )
 9602 format('>> Terminating analysis at this point',//)
 9900 format('>> incremental displacements: ',
     & 1000(/4x,i8, f12.8) )
c
c
       contains
c      ========
c
c
c     ****************************************************************
c     *                                                              *
c     *   set up equation sparsity for local assembly: symmetric     *
c     *                                                              *
c     *                       written by  : rhd                      *
c     *                   last modified : 11/19/2019 rhd
c     *                                                              *
c     ****************************************************************
c
      subroutine ds_setup_sparsity_local( ireturn, using_mpi,
     &                                    num_ranks )
      implicit none
c
      integer :: ireturn, num_ranks  ! local to this routine
      logical :: using_mpi
c
c                    locals
c
      integer :: i, j, k, l, srow, start_srow
      integer, allocatable :: start_kindex_locs(:), edest(:,:,:),
     &                        scol_flags(:,:), scol_lists(:,:)
      character(len=200) mkl_string
      integer :: mkl_num_thrds, next_space, count_previous, count_now,
     &           nrow_lists, safety_factor
      logical, parameter :: local_debug_2 = .false.
c
c              1. generate equation numbers for all unconstrained
c                 nodal dof. dof_eqn_map(i) gives equation
c                 number for structure dof number. eqn_node_map
c                 gives the node number for an equation number.
c                 set the actual number of equations we will solve
c                 (neqns). if this turns out = 0, all nodes in
c                 model are fully constrained and we just leave.
c                 the solution is known from constraints.
c
c
      if( local_debug ) write(out,*) '.. @ (1)  '
      call thyme( 21, 1 )
      if( cpu_stats .and. show_details ) then
          call mkl_get_version_string( mkl_string )
          mkl_num_thrds = solver_threads
          if( using_mpi ) then
             write(out,9300) mkl_num_thrds, num_ranks,
     &                       mkl_string(38:45),
     &                       mkl_string(61:68), wcputime(1)
           else
             write(out,9400) mkl_num_thrds, mkl_string(38:45),
     &                       mkl_string(61:68), wcputime(1)
          end if
      end if
      if( .not. allocated ( dof_eqn_map ) ) then
          allocate( dof_eqn_map(num_struct_dof) )
          allocate( eqn_node_map(num_struct_dof) )
      end if
      call dof_map( dof_eqn_map, cstmap, nonode, num_enode_dof,
     &              eqn_node_map, neqns )
      if( cpu_stats .and. show_details ) write(out,9409) wcputime(1)
c
      ireturn = 1
      if( neqns .le. 0 ) return
c
c              2. decide if this set of equations has a different
c                 sparsity from the previous time we came thru here.
c                 right now, this is based on the actual number
c                 of equations. in the future, this test may need
c                 to actual verify the same sparsity structure, e.g.,
c                 the same numberof equations and the same
c                 number of off diagonal terms.
c
      if( neqns .ne. old_neqns )  new_size = .true.
      if( local_debug ) write(out,*) '....  neqns, old_neqns: ',
     &                             neqns, old_neqns
      old_neqns  = neqns
      ireturn = 2
      if( .not. new_size ) return
c
c              3. compute the number of non-zero terms in the
c                 upper-triangle of the assembled structure stiffness.
c                 this is done by simulating a row-by-row
c                 sparse format assembly (rather than an element-by-
c                 element assembly). this way we never need to allocate
c                 memory for the profile of the upper-triangle.
c                 we save the count of the number terms in the
c                 "assembled" profile - the count can become
c                 changed later in this routine if, for example, there
c                 are mpc equations which modify the starting
c                 "assembled" stiffness.
c
c                 we also get the number of non-zero terms on each row
c                 to right of diagonal ( start_kindex_locs )
c
      safety_factor = 10
      allocate( start_kindex_locs(neqns) )
      allocate( scol_flags(neqns,num_threads) )
      allocate( edest(mxedof,mxconn,num_threads) )
      edest = 0
      nrow_lists = mxconn * mxedof * safety_factor
      allocate( scol_lists(nrow_lists,num_threads) )
c
      call count_profile_symmetric( neqns, ncoeff, num_threads,
     &                              eqn_node_map, iprops, dof_eqn_map,
     &                              start_kindex_locs, edest,
     &                              scol_flags, scol_lists,
     &                              nrow_lists  )
c
      num_terms  = neqns + ncoeff
      old_ncoeff = ncoeff
      ncoeff_from_assembled_profile = ncoeff
      call thyme( 21, 2 )
      call thyme( 22, 1 )
c
      if( cpu_stats .and. show_details ) then
          write(out,9410) wcputime(1)
          write(out,9411) neqns
          write(out,9412) num_terms
      end if
c
c              4. build sparsity of equilibrium equations
c                 as defined by the k_indexes, k_ptrs
c                 vectors. see diagram at end for example.
c                 The work is done is sparse format so that
c                 we only need a workspace vector of neqns
c                 in length. Refer to comments inside called
c                 routine.
c
c                 start_kindex_locs from count.. routine has k_ptrs.
c                 convert to the starting locating in k_indexes for
c                 storage of column indexes for each equation.
c
      start_kindex_locs(neqns) = 0
c
      if( local_debug_2 )  then
          write(out,*) '...  start_kindex_locs ...'
          do i = 1, neqns
             write(out,*) '        ',i,start_kindex_locs(i)
          end do
      end if
c
      select case( neqns )
        case( 1 )
        case( 2 )
          start_kindex_locs(1) = 1
          if( ncoeff_from_assembled_profile == 0 ) ! 2x2 no (1,2) term
     &           start_kindex_locs(1) = 0
        case( 3: )
          start_srow = 0
          do srow = 1, neqns-1
            if( start_kindex_locs(srow) > 0 ) then
               start_srow = srow
               exit
            end if
          end do
          if( start_srow == 0 ) then
               write(out,9210)
               call die_abort
          end if
          count_previous = start_kindex_locs(start_srow)
          start_kindex_locs(start_srow) = 1
          do i = start_srow+1, neqns-1
           count_now = start_kindex_locs(i)
           start_kindex_locs(i) = start_kindex_locs(i-1) +
     &                            count_previous
           count_previous = count_now
          end do
          if( local_debug_2 ) then
              write(out,*) '... start_k_indexes_locs after update'
              do i = 1, neqns
                 write(out,*) '        ',i,start_kindex_locs(i)
              end do
          end if
      end select
c
      if( cpu_stats .and. show_details ) write(out,9419) wcputime(1)
      if( allocated( save_k_indexes ) ) deallocate( save_k_indexes )
      if( allocated( save_k_ptrs ) )    deallocate( save_k_ptrs )
      allocate( save_k_indexes(ncoeff_from_assembled_profile) )
      allocate( save_k_ptrs(neqns) )
      if( local_debug ) write(out,*)
     &                 ' @ 2.2 ... saving K-indexes, ptrs'
c
      call build_col_sparse_symm(
     &      neqns, num_threads, eqn_node_map, dof_eqn_map,
     &      save_k_indexes, save_k_ ptrs,
     &      start_kindex_locs, edest, scol_lists, nrow_lists, ncoeff )
c
      deallocate( start_kindex_locs, scol_flags, edest, scol_lists )

      if( cpu_stats .and. show_details ) write(out,9420) wcputime(1)
      call thyme( 22, 2 )
c
      ireturn = 3
      return
c
 9210 format('>> FATAL ERROR: computations of dense storage data',
     & /,    '                failed in equation solver.',
     & /,    '                inconsistencies in data structure',
     & /,    '                job terminated' )
 9300  format (
     &  10x, '>> solver wall time statistics (secs):'
     & /,15x,'number of MKL threads used      ',i10,
     & /,15x,'number of MPI ranks:            ',i10,
     & /,15x,'MKL version, build: ',5x,a8,1x,a8,
     & /,15x,'starting work                 @ ',f10.2 )
 9400  format (
     &  10x, '>> solver wall time statistics (secs):'
     & /,15x,'number of MKL threads used      ',i10,
     & /,15x,'MKL version, build: ',5x,a8,1x,a8,
     & /,15x,'starting work                 @ ',f10.2 )
 9409  format(
     &  15x, 'finished dof setup            @ ',f10.2 )
 9410  format(
     &  15x, 'non-zero term counts finished @ ',f10.2 )
 9411  format(
     &  15x, 'number of equations             ',i10)
 9412  format(
     &  15x, 'non-zero terms in profile       ',i10)
 9419  format(
     &  15x, 'building sparse ptrs, indexes @ ',f10.2 )
 9420  format(
     &  15x, 'k_ptrs, k_indexes done        @ ',f10.2 )
c
      end subroutine ds_setup_sparsity_local

c     ****************************************************************
c     *                                                              *
c     *     set up equation sparsity for distributed assembly        *
c     *                                                              *
c     *                       written by  : rhd                      *
c     *                   last modified : 04/13/2015                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine ds_setup_sparsity_distributed( ireturn )
      implicit none
c
      integer :: ireturn  ! local to this routine
c
c           If we have the same sparsity pattern we can skip
c           all of this and simply assemble the coefficients
c
      if( .not. allocated ( dof_eqn_map ) ) then
          allocate( dof_eqn_map(num_struct_dof) )
          allocate( eqn_node_map(num_struct_dof) )
      end if
      call dof_map( dof_eqn_map, cstmap, nonode, num_enode_dof,
     &              eqn_node_map, neqns )
c
      ireturn = 1
      if( neqns .le. 0 ) return
c
      if( neqns .ne. old_neqns ) new_size = .true.
      old_neqns = neqns
      deallocate(dof_eqn_map)
      deallocate(eqn_node_map)
c
      ireturn = 2
      if( .not. new_size ) return
c
c           Determine the sparsity structure of your blocks of elements
c
      call wmpi_alert_slaves( 41 )
      call determine_local_sparsity
c
c           Determine the initial map for assembling the sparse
c           structure this should just be a simple heuristic to keep
c           a good amount of locality
c
      call wmpi_alert_slaves( 42 )
      call determine_initial_map
c
c           Assemble the sparse structure for the above map, now we know
c           global nnz and full sparsity.  We can make a parmetis graph
c           to balance.
c
      call wmpi_alert_slaves( 43 )
      call assemble_sparsity
c
c           Determine a better mapping for load balance, pass rows as
c           necessary and make a note of the new equation numbers
c
      call wmpi_alert_slaves( 44 )
      call determine_ordering
c
c           Move the nnz structure around to match the final map
c
      call wmpi_alert_slaves( 45 )
      call move_sparsity
c
c           all done
c
      ireturn = 3
c
      return
c
      end subroutine ds_setup_sparsity_distributed

c     ****************************************************************
c     *                                                              *
c     *                        ds_root_assembly                      *
c     *                                                              *
c     *                       written by  : rhd                      *
c     *                   last modified : 09/14/2015                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine ds_root_assembly
      implicit none
c
      logical :: hypre_solver, rebuild_mpcs
      double precision :: px, py, pz
c
c              1. sparsity of equilibrium equations already defined on
c                 this or previous solution. neqns = 0 means all dof in
c                 model are constrained. just leave.
c
      if( local_debug ) write(out,*) ' @ 3'
      if( neqns .le. 0 ) return
c
c              2. make a working copy of the k_ptrs and k_indexes from
c                 saved values for use by the sparse solvers
c                 (some sprase solvers may alter them and/or
c                 the values/sizes are changed during insertion of
c                 mpcs). allocate space for
c                 the diagonal terms, the off-diagonal, non-zero
c                 coefficients and the effective load vector. initialize
c                 terms as needed.
c
c                 The dimension of the arrays are inflated in order to
c                 allow using different sparse solvers. The format of
c                 these vectors in the assembly process follow
c                 arrays in the now old old NASA-VSS sparse solver.
c
c                 Other solvers
c                 will have to perform mapping/transformation of these
c                 arrays to their native format, e.g if the diagonal
c                 terms are snot stored in a separate vector.
c                 Therefore, the sizes
c                 of these arrays are larger than the needed sizes
c                 in the NASA-VSS solver.
c
c                 NASA-VSS solver uses:
c
c                       k_ptrs(neqns)
c                       k_indexes(ncoeff)
c                       k_diag(neqns)
c                       k_coeffs(ncoeff)
c                       p_vec(neqns)
c                  u_vec(neqns)
c
c                 The hypre solver needs a non-symmetric CSR
c                 format matrix.  The hypre functions take care of
c                 conversion, but we need to overallocate the
c                 arrays in order to fit all the
c                 data.  Note for non-symmetric CSR we need:
c
c                       k_ptrs(neqns+1)
c                       k_indexes(2*ncoeff+neqns)
c                       k_coeffs(2*ncoeff+neqns)
c                       p_vec(neqns)
c                       u_vec(neqns)
c
c                 k_diag is not used in the CSR format, but we'll
c                 need it for assembly.
c
      hypre_solver =  solver_flag .eq. 9
      call thyme( 18, 1 )
      if( hypre_solver .or. asymmetric_assembly  ) then
            allocate( k_ptrs(neqns + 1 ),
     &                k_indexes(2*ncoeff + neqns),
     &                k_coeffs(2*ncoeff + neqns) )
      else
            allocate( k_ptrs(neqns + 1),
     &                k_indexes(ncoeff + neqns),
     &                k_coeffs(ncoeff + neqns) )
      end if
      allocate( k_diag(neqns), p_vec(neqns), u_vec(neqns) )
      if( local_debug ) write(out,*) '  @ 3.1'
      do i = 1, neqns
        k_ptrs(i) = save_k_ptrs(i)
        k_diag(i) = zero
        p_vec(i)  = zero
        u_vec(i)  = zero
      end do
      if( local_debug ) then
         write(out,*)  '  @ 3.2, ncoeff: ',ncoeff
         write(out,*)  '         ncoeff_from_assembled_profile: ',
     &                        ncoeff_from_assembled_profile
         write(out,*)  '         size k_indexes: ', size(k_indexes)
      end if
      k_indexes(1:ncoeff_from_assembled_profile) =
     &            save_k_indexes(1:ncoeff_from_assembled_profile)
      if( local_debug ) write(out,*) ' @ 4'
c
c              3. Here we branch for full (asymmetric) or
c                 upper-triangle symmetric assembly.
c
c                 the code in assem_by_row runs parallel with threads
c                 where a thread builds an entire equation row
c
      if( asymmetric_assembly ) then
c
c             3a. The only difference is we convert the UT
c                 VSS ptrs and indexes to a full, structurally
c                 symmetric CSR format and use a modified
c                 assem_by_row which picks up the LT terms.
c                 rows are assembled in parallel
c
        nnz = 0  ! total number of non-zero terms in assembled eqns.
        call convert_vss_csr( neqns, ncoeff, save_k_ptrs(1),
     &      save_k_indexes(1), nnz, k_ptrs(1), k_indexes(1) )
        if( cpu_stats .and. show_details ) write(out,9999) wcputime(1)
        k_coeffs = zero
        call assem_by_row_asymmetric( neqns, num_threads, eqn_node_map,
     &                                dof_eqn_map, k_ptrs, k_indexes,
     &                                iprops, k_coeffs )

      else
c             3b. assemble symmetric  equilibrium equations
c                 in sparse format using a row-by-row algorithm
c                 rather than a conventional element-by-element.
c                 rows are assembled in parallel.
c
        k_coeffs = zero
        call assem_by_row( neqns, num_threads, eqn_node_map,
     &                     dof_eqn_map, k_diag, k_coeffs,
     &                     k_indexes, k_ptrs, iprops,
     &                     dcp, noelem )
c
      end if ! for asymmetric/symmetric assembly
c
      if( cpu_stats .and. show_details ) write(out,9470) wcputime(1)
      if( local_debug ) write(out,*) ' @ 5'
c
c              4. set the rhs of the equations. p_vec has the
c                 neqns terms for unconstrained dof from res.
c
      if( local_debug3 ) then
         write(out,9500)
         i_offset = 0
         do i = 1, nonode
           px = res(i_offset+1)
           py = res(i_offset+2)
           pz = res(i_offset+3)
           i_offset = i_offset + 3
           write(out,9510) i, px, py, pz
         end do
      end if
c
      do i = 1, nodof
        eqn_num = dof_eqn_map(i) ! dof # -> eqn # map
        if( eqn_num .ne. 0 ) p_vec(eqn_num) = res(i)
      end do
c
      call thyme( 18, 2 )
      if( local_debug ) write(out,*) ' @ 6'
      if( cpu_stats .and. show_details ) write(out,9480) wcputime(1)
c
c              5. if any multi-point constraint equations exist,
c                 modify the stiffness matrix and rhs accordingly.
c                 if the # eqns has not changed, a much faster set
c                 of routines runs, taking advantage of info learned
c                 during the previous solve. note that the mpc
c                 insertion process for a new set of equations
c                 (new_size is .true.) resets the value of ncoeff !
c                 Thats why we save the "assembled" profile count
c                 in step 3 of the process above to construct the
c                 equation sparsity.
c
c                 We also need to recreate MPCs if they've
c                 been modified! Added a flag for this case
c
c                 ==>> MPCs/tied cons are supported only for symmetric
c                      assembly
c
c                 The re-insert code has a bug. minimize effects
c                 by rebuilding MPC data structures on 1st iteration
c                 of each step.
c
      if( local_debug2 ) then
          write(out,*) '... before MPC enforcement ...'
          write(out,*) '       neqns, ncoeff: ', neqns, ncoeff
          write(out,*) '       .... k_ptrs ....'
          write(out,9010) (i, k_ptrs(i), i = 1, neqns+1)
          k_ptrs(neqns+1) = 0
          write(out,*) '       .... k_indexes ....'
          write(out,*) '      sizeof k_indexes: ', sizeof(k_indexes)
          write(out,9010) (i, k_indexes(i), i = 1, ncoeff+neqns)
      end if

      if( tied_con_mpcs_constructed .or. mpcs_exist ) then
          rebuild_mpcs = now_iteration .eq. 1
          if( new_size .or. modified_mpcs .or. rebuild_mpcs ) then
            if( local_debug ) write(out,*) '.... @ 6.1'
            call mpc_insert_terms(neqns, k_ptrs, k_diag, dstmap,
     &                            dof_eqn_map)
            if( local_debug ) write(out,*) '.... @ 6.2'
            call mpc_modify_stiffness(neqns, k_diag, p_vec)
            if( local_debug ) write(out,*) '.... @ 6.3'
            call mpc_remove_dep_eqns(neqns, k_ptrs)
            if( local_debug ) write(out,*) '.... @ 6.4'
            modified_mpcs = .false.
         else
            if( local_debug ) write(out,*) '.... @ 6.5'
            call mpc_reinsert_terms(neqns, k_ptrs, dstmap, dof_eqn_map)
            if( local_debug ) write(out,*) '.... @ 6.6'
            call mpc_modify_stiffness(neqns, k_diag, p_vec)
            if( local_debug ) write(out,*) '.... @ 6.7'
            call mpc_remove_dep_eqns(neqns, k_ptrs)
            if( local_debug ) write(out,*) '.... @ 6.8'
         end if
         if( cpu_stats .and. show_details ) then
            write(out,9490) wcputime(1)
            write(out,9412) ncoeff+neqns
         end if
         if( local_debug ) then
            write(out,*) "... pvec after MPC mods ..."
            do k = 1, neqns
               write(out,9300) k, p_vec(k)
            end do
         end if

      end if ! on tied_con_mpcs_constructed .or. mpcs_exist
c
      if( local_debug2 ) then
          write(out,*) '... after MPC enforcement ...'
          write(out,*) '       neqns, ncoeff: ', neqns, ncoeff
          write(out,*) '       .... k_ptrs ....'
          write(out,9010) (i, k_ptrs(i), i = 1, neqns+1)
          write(out,*) '       .... k_indexes ....'
          write(out,9010) (i, k_indexes(i), i = 1, ncoeff+neqns)
      end if
c
      if( local_debug ) write(out,*) ' @ 7'

c
c              6. write equations ready to solve in various formats
c                 as requested by user. makes our equations available
c                 to drive solvers w/o dragging warp3d along
c
      if( sparse_stiff_output ) then
          call save_solver_sparse( neqns, ncoeff,
     &                           k_diag, p_vec, k_coeffs, k_ptrs,
     &                           k_indexes, dstmap, dof_eqn_map,
     &                           eqn_node_map, num_struct_dof,
     &                           nonode, sparse_stiff_output,
     &                           sparse_stiff_binary,
     &                           sparse_stiff_file_name, out )
          if( sparse_research )
     &         call store_solver( sparse_stiff_binary,
     &                            sparse_stiff_file_name,
     &                            neqns, num_struct_dof,
     &                            dof_eqn_map, eqn_node_map )
      end if
c
c              7. deallocate all element stiffness matrices.
c                 6/23/2017 - we now run w/o deleteing all those
c                 blocks on each global Newton iteration.
c
c      call estiff_allocate( 5 )
      if( local_debug ) write(*,*) ' @ 8'
c
      return
c
 9010 format(6x,12i8)
 9300 format(10x,i5,3(1x,f10.6))
 9412  format(
     &  15x, 'non-zero terms in profile       ',i10)
 9470  format(
     &  15x, 'sparse [k] assembly done      @ ',f10.2 )
 9480  format(
     &  15x, 'load vector assembly done     @ ',f10.2 )
 9490  format(
     &  15x, 'multi-point constraints added @ ',f10.2 )
 9500 format( 2x,"... Force vector for equation solving: ",
     &   /,   2x,"   node       /----   px, py, pz    ----/")
 9510 format( 2x,i7, 3e14.6 )
 9999 format(
     &  15x, 'upper triangle to full done   @ ',f10.2)
c
      end subroutine ds_root_assembly

c     ****************************************************************
c     *                                                              *
c     *                     ds_distributed_assembly                  *
c     *                                                              *
c     *                       written by  : rhd                      *
c     *                   last modified : 04/14/2015                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine ds_distributed_assembly
      implicit none
c
c
c           The idea here is to move the sparsity to the correct processors
c           and then assemble the actual data.
c           The first part (the sparsity) would be done above and would save
c           memory, but because we might spawn processors must be done here.
c           Move the nnz structure around to match the final map
c
      call wmpi_alert_slaves( 46 )
      call assemble_coefs( new_size )
c
c           and get and distribute the load vector
c
      call wmpi_alert_slaves( 47 )
      call assem_load_vec
      call wmpi_alert_slaves( 48 )
      call dist_final_setup
c
c           We need to allocate the solution vector on root
c
      neqns = curr_neqns()   !   a function in mpi code
      allocate( u_vec(neqns) )
c
c           It turns out that we need this
c
      if( .not. allocated ( dof_eqn_map ) ) then
          allocate( dof_eqn_map(num_struct_dof) )
          allocate( eqn_node_map(num_struct_dof) )
      end if
      call dof_map( dof_eqn_map, cstmap, nonode, num_enode_dof,
     &              eqn_node_map, neqns )
c
      return
c
c           Checksum matches.  We have the right distributed equations
c
      end subroutine ds_distributed_assembly

c     ****************************************************************
c     *                                                              *
c     *                     ds_drive_hypre                           *
c     *                                                              *
c     *                       written by  : rhd                      *
c     *                   last modified : 04/14/2015                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine ds_drive_hypre
      implicit none

       if( asymmetric_assembly ) then
          write(out,*) "Note to Mark.  This should work, but check"
        end if
c
        call wmpi_alert_slaves( 40 )
c
c           run the hyper solver on MPI ranks
c
        call iterative_sparse_hypre( u_vec, out, error_code )
c
c           The returned error_code means one of two things:
c
c           1. we've run into the condition "hypre can't solve the
c              increment before an adaptive reset is necessary", or
c
c           2. "these equations are too ill-conditioned for hypre
c              to solve."
c
c           If the adaptive Newton solution flag is not set
c           assume (2) and just exit.
c
c           Else hypre will keep count of how many adaptive
c           calls in a row it has triggered.  If this is the
c           first hypre-triggered adaptive step, then deallocate
c           the structs and return to mnalgr, requesting a usual Newton
c           adaptive sub-step.  If this is the second (-in a row-),
c           then quit out, as we're likely not going to recover.
c
      if( error_code .eq. 1 ) then
         if( adaptive_flag ) then
           if( precond_fail_count .ge. 2 ) then
              write(out,9600)
              write(out,9602)
              call die_gracefully
           else
              if( local_debug ) write (out,*) '... @ 8.1'
              write(out,9605)
              hyp_trigger_step = .true.
              if( .not. (parallel_assembly_used) ) then
                deallocate( k_diag )
                deallocate( k_coeffs )
                deallocate( k_ptrs )
                deallocate( k_indexes )
                deallocate( p_vec )
              end if
              deallocate( u_vec )
              return
           end if
         else ! user does not allow adaptive Newton.
               write(out,9607)
               write(out,9602)
               call die_gracefully
         end if
      end if
c
c           this means we solved the equations and do not need an
c           adaptive call, so set the flag accordingly.
c
      hyp_trigger_step = .false.
c
      return
c
 9600 format('>> ERROR: hypre has requested 2 adaptive steps',
     &          ' in a row due to non-convergence.')
 9602 format('>> Terminating analysis at this point',//)
 9605 format('>> Note: hypre has requested adaptive step reduction'/)
 9607 format('>> ERROR: hypre failed to converge. Try adjusting load',
     &  /    '          step size on hypre parameters.'//)
c
      end subroutine ds_drive_hypre

      end subroutine drive_assemble_solve

c ------------------------------------------------------------------------
c                        NASA - VSS
c
c               sparse matrix storage format
c
c   example:
c
c                1    2    3    4     5     6
c
c          1  | 100   1    2                5  |  | x1 |     | 201 |
c          2  |     200    6    7           9  |  | x2 |     | 202 |
c          3  |          300   10    11    12  |  | x3 |     | 203 |
c      a = 4  |                400   13    14  |  | x4 |  =  | 204 |
c          5  |                     500    15  |  | x5 |     | 205 |
c          6  |                           600  |  | x6 |     | 206 |
c
c     number of equations    = 6
c
c     number of coefficients = 12
c
c
c     k_ptrs    = { 3, 3, 3, 2, 1, 0}
c
c     k_indesxs = { 2, 3, 6,  3, 4, 6,  4, 5, 6,  5, 6,  6}
c
c     k_coefs   = { 1, 2, 5,  6, 7, 9, 10,11,12, 13,14, 15}
c
c     k_diag    = { 100, 200, 300, 400, 500, 600}
c
c     k_rhs     = { 201, 202, 203, 204, 205, 206}
c
c ------------------------------------------------------------------------
c     ****************************************************************
c     *                                                              *
c     *           subroutines save_solver_sparse                     *
c     *           write the equations in sparse format onto a        *
c     *           stream, unformatted file or a formatted file       *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 02/15/2014                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine save_solver_sparse( neqns, ncoeff, k_diag, p_vec,
     &    k_coeffs, k_ptrs, k_indexes, dstmap, dof_eqn_map,
     &    eqn_node_map, num_struct_dof, nonode, sparse_stiff_output,
     &    sparse_stiff_binary, sparse_stiff_file_name, iout )
c
      implicit none
      integer  neqns, ncoeff, k_ptrs(*), k_indexes(*), dstmap(*),
     &         dof_eqn_map(*), eqn_node_map(*), num_struct_dof, nonode
      double precision
     &  k_diag(*), p_vec(*), k_coeffs(*)
      logical sparse_stiff_output, sparse_stiff_binary, connected
      character(len=*) :: sparse_stiff_file_name
      integer fileno, open_result, iout, i
c
      sparse_stiff_output = .false.
c
c                  find an available unit number to use
c
      do fileno = 11, 99
        inquire(unit=fileno, opened=  connected )
        if ( .not. connected ) go to 100
      end do
      write(iout,9000)
      return
c
 100  continue
      if ( sparse_stiff_binary ) then
        open(unit=fileno, file=sparse_stiff_file_name,
     &       status='replace', access='stream',form='unformatted',
     &       iostat=open_result  )
      else
        open(unit=fileno,file=sparse_stiff_file_name,
     &       status='replace', access='sequential', form='formatted',
     &       iostat=open_result  )
      end if
      if ( open_result .ne. 0 ) then
         write(iout,9010)
         return
      end if
c
      if ( sparse_stiff_binary ) then
          write(iout,9020)
          write(fileno) neqns, ncoeff
          write(iout,9022)
          write(fileno) (k_diag(i), i=1,neqns)
          write(iout,9024)
          write(fileno) (p_vec(i), i=1,neqns)
          write(iout,9026)
          write(fileno) (k_coeffs(i), i=1,ncoeff)
          write(iout,9028)
          write(fileno) (k_ptrs(i), i=1,neqns)
          write(iout,9030)
          write(fileno) (k_indexes(i), i=1,ncoeff)
          write(iout,9032)
          write(fileno) num_struct_dof, nonode
          write(fileno) dstmap(1:nonode)
          write(fileno) dof_eqn_map(1:num_struct_dof)
          write(fileno) eqn_node_map(1:neqns)
          close( fileno,status='keep' )
          write(iout,9034)
          return
      end if
c
      write(iout,9020)
      write(fileno,9500) neqns, ncoeff
      write(iout,9022)
      write(fileno,9510) (k_diag(i), i=1,neqns)
      write(iout,9024)
      write(fileno,9510) (p_vec(i), i=1,neqns)
      write(iout,9026)
      write(fileno,9510) (k_coeffs(i), i=1,ncoeff)
      write(iout,9028)
      write(fileno,9520) (k_ptrs(i), i=1,neqns)
      write(iout,9030)
      write(fileno,9520) (k_indexes(i), i=1,ncoeff)
      write(fileno,9500) num_struct_dof, nonode
      write(fileno,9520) dstmap(1:nonode)
      write(fileno,9520) dof_eqn_map(1:num_struct_dof)
      write(fileno,9520) eqn_node_map(1:num_struct_dof)
      write(iout,9032)
      close( fileno,status='keep' )
      write(iout,9034)
      return
c
 9000 format('>> WARNING: could not find a unit number to write',
     & /     '            sparse matrix. action skipped...',/)
 9010 format('>> WARNING: could not open file to save sparse [K]',
     & /     '            action skipped...',/)
 9020 format(15x,'writing sparse stiffness in binary format...')
 9022 format(17x,'> matrix sizes written...')
 9024 format(17x,'> diagonal terms written...')
 9026 format(17x,'> rhs written...')
 9028 format(17x,'> non-zero upper triangle terms written...')
 9030 format(17x,'> row pointers written...')
 9032 format(17x,'> non-zero column numbers written...')
 9034 format(17x,'> file closed...')
 9500 format(2i10)
 9510 format(4e20.12)
 9520 format(7i12)
c
      end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine store_solver                 *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 1/15/02 rhd                *
c     *                                                              *
c     *     write a binary file of key structure data after assembly *
c     *     and before solve by sparse solver. this file is used to  *
c     *     support external programs for solver research.           *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine store_solver( do_binary, file_name, neqns,
     &                         num_struct_dof, dof_eqn_map,
     &                         eqn_node_map )
      use global_data ! old common.main
      use main_data
      implicit integer (a-z)
      intrinsic size
      logical do_binary
      character(len=*) :: file_name
      dimension dof_eqn_map(*), eqn_node_map(*)
c
c               parameter and local variable declarations
c
      allocatable ele_info(:)
      real dumr, rzero
      double precision
     &     dumd
      character(len=1) :: dums
      data check_data_key, rzero / 2147483647, 0.0 /
c
      if ( .not. do_binary ) then
         write(out,*) '>>> WARNING: solver research data cannot'
         write(out,*) '             be written to an ascii file'
         write(out,*) '             request ignored....'
         return
      end if
c
      prec_fact = 2
      fileno = 11
        open(unit=fileno, file=file_name,
     &       status='old', access='sequential', form='unformatted',
     &       iostat=open_result, position='append'  )
      if ( open_result .ne. 0 ) then
         write(out,9070)
         return
      end if
      write(out,9050)
c
c                       write out sizes and integer scalars
c
      write(fileno) noelem, nonode, nodof, csthed, inctop, nelblk
      write(fileno) check_data_key
      write(out,9000)
c
c                       write out integer vectors
c
      call wrtbk( fileno, invdst, nodof )
      call wrtbk( fileno, incmap, noelem )
      call wrtbk( fileno, dstmap, nonode )
      call wrtbk( fileno, cstmap, nodof )
      call wrtbk( fileno, incid,  inctop )
      allocate( ele_info(noelem) )
      do i = 1, noelem
        ele_info(i) = iprops(2,i)
      end do
      call wrtbk( fileno, ele_info, noelem )
      deallocate( ele_info )
      write (fileno) check_data_key
c
c                       write more integer vectors, arrays
c
      call wrt2d( fileno, elblks(0,1) , 4, 4, nelblk  )
      call wrtbk( fileno, cp, mxedof )
      call wrtbk( fileno, dcp, mxedof )
      call wrt2d( fileno, icp, mxutsz, mxutsz, 2  )
      write (fileno) check_data_key
      write(out,9010)
c
c                       constraint data (double)
c
      call wrtbk( fileno, cnstrn, prec_fact*nodof )
      write (fileno) check_data_key
      write(out,9040)
c
c                       inverse incidences
c
      call store_solver2( fileno )
      write (fileno) check_data_key
c
c                       element stiffness blocks and element
c                       destination vectors
c
      call store_solver3( fileno )
      write(out,9080)
c
c                       number of actual equations to solve, the
c                       structure dof -> equation no. map vector, and
c                       equation no. -> structure node map vector
c
      write(fileno) neqns, num_struct_dof
      write(fileno) dof_eqn_map(1:num_struct_dof)
      write(fileno) eqn_node_map(1:num_struct_dof)
      write(fileno) check_data_key
c
c                       close the file
c
      close( fileno, status='keep' )
      write(out,9060)
c
c
c
 1000 format ( 3x, e16.6 )
 9000 format(17x,'> scalars written...')
 9010 format(17x,'> integer arrays written...')
 9040 format(17x,'> double precision arrays written...')
 9050 format(15x,'appending solver research data to binary file')
 9060 format(15x,'solver research file written and closed...')
 9070 format('>> WARNING: could not open file to save sparse [K]',
     & /     '            action skipped...',/)
 9080 format(17x,'> element stiffness and destinations written...')

      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine store_solver_2                    *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 01/15/02                   *
c     *                                                              *
c     *     write data into restart file for the requested data      *
c     *     structures which additional complexity                   *
c     *                                                              *
c     ****************************************************************
c
      subroutine store_solver2( fileno )
      use global_data ! old common.main
c
      use main_data, only : inverse_incidences, inverse_dof_map
c
      implicit integer (a-z)
      data check_data_key / 2147483647 /
c
      write(fileno) ( inverse_incidences(stnd)%element_count, stnd = 1,
     &                nonode )
      do stnd = 1, nonode
        ecount = inverse_incidences(stnd)%element_count
        write(fileno) ( inverse_incidences(stnd)%element_list(i),
     &                  i = 1, ecount )
      end do
      write(out,9000)
c
      do stnd = 1, nonode
        ecount = inverse_incidences(stnd)%element_count
        write(fileno) (( inverse_dof_map(stnd)%edof_table(i,j),
     &                  i = 1, ecount ), j = 1, 3 )
      end do
      write(out,9100)
c
      return
 9000 format(17x,'> inverse_incidences written...')
 9100 format(17x,'> inverse_dof_maps written...')
      end
c     ****************************************************************
c     *                                                              *
c     *                 subroutine store_solver_3                    *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 01/15/02                   *
c     *                                                              *
c     *     write data into restart file for the requested data      *
c     *     structures which additional complexity                   *
c     *                                                              *
c     ****************************************************************
c
      subroutine store_solver3( fileno )
      use global_data ! old common.main
c
      use elem_block_data, only : edest_blocks, estiff_blocks
c
      implicit integer (a-z)
      data check_data_key / 2147483647 /
c
      do blk = 1, nelblk
        felem         = elblks(1,blk)
        span          = elblks(0,blk)
        nnode         = iprops(2,felem)
        num_enode_dof = iprops(4,felem)
        totdof        = nnode * num_enode_dof
        write(fileno) nnode, totdof, span
        write(fileno)  edest_blocks(blk)%ptr
        write(fileno)  estiff_blocks(blk)%ptr
        write(fileno) check_data_key
      end do
      return
      end

c     ****************************************************************
c     *                                                              *
c     *                 subroutine store_csr_formatted               *
c     *                                                              *
c     *                       written by : mcm                       *
c     *                                                              *
c     *                   last modified : 12/18/2016                 *
c     *                                                              *
c     *     write assembled equations IN CSR in the following format:*
c     *     n, nnz  (2I8)                                            *
c     *     value(i), col_index(i) (E25.16, 2X, I8)                  *
c     *     row_ptrs (I8)                                            *
c     *     RHS (E25.16)                                             *
c     *                                                              *
c     ****************************************************************
c
      subroutine store_csr_formatted( filename, n, nnz, row_ptrs,
     &                        col_indexes, values, rhs, diagonal )
      implicit none

      integer :: n, nnz, row_ptrs(*), col_indexes(*)
      double precision ::  diagonal(*), rhs(*), values(*)
      character(len=*) :: filename

      integer :: i, fileno
      logical :: connected
c
c          Find a file number
c
      do fileno = 11, 99
        inquire(unit=fileno, opened=connected )
        if ( .not. connected ) go to 100
      end do
      write(*,*) "Couldn't find a file number to write to, skipping"
      return
c
 100  continue
c
c           Basically the plan is to convert to CSR (with existing fn)
c           and write out in
c           the format specified above.
c
      call pardiso_symmetric_map( n, nnz, diagonal, rhs,values,
     &                   row_ptrs, col_indexes )
      nnz = nnz + n
      open (unit = fileno, file=filename)
      write(*,*) "Writing"
      write(fileno,5) n,nnz
5     format(2I10)
      do i=1, nnz
            write(fileno,10) values(i), col_indexes(i)
      end do
10    format(E25.16,2X,I10)
      do i = 1, n+1
            write(fileno,15) row_ptrs(i)
      end do
15    format(I10)
      do i = 1, n
            write(fileno,20) rhs(i)
      end do
20    format(E25.16)
      close(fileno)
      write (*,*) "Done writing"
      return
      end
c
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine convert_vss_csr                   *
c     *                                                              *
c     *                       written by : mcm                       *
c     *                                                              *
c     *                   last modified : 12/10/13                   *
c     *                                                              *
c     *     Convert the upper triangle vss sparsity pattern to       *
c     *     a full CSR format.                                       *
c     *                                                              *
c     ****************************************************************
c
      subroutine convert_vss_csr( neqns, ncoeff, vss_ptrs,
     &                            vss_indexes, nnz, ptrs, indexes )
      implicit none

      integer, intent(in) :: neqns, ncoeff, vss_ptrs(neqns+1),
     &                        vss_indexes(neqns+ncoeff)
      integer, intent(out) :: nnz, ptrs(neqns+1),
     &                        indexes(2*ncoeff+neqns)
      integer :: i, j, s, o, k, c, r, work(neqns)
c
      nnz = neqns + 2*ncoeff
c
c           Copy the number of terms per row
c
      ptrs = 0
      ptrs(1:neqns) = vss_ptrs(1:neqns)
c
c           Add the diagonal terms
c
      ptrs(1:neqns) = ptrs(1:neqns) + 1
c
c           Add the number of terms per column, store for later use
c
      work = 0
      do i = 1, ncoeff
        ptrs(vss_indexes(i)) = ptrs(vss_indexes(i)) + 1
        work(vss_indexes(i)) = work(vss_indexes(i)) + 1
      end do
c
c           Convert to offsets
c
      do i = 2, neqns
        ptrs(i) = ptrs(i-1) + ptrs(i)
      end do
      do i = neqns+1, 2, -1
        ptrs(i) = ptrs(i-1) + 1
      end do
      ptrs(1) = 1
c
c           Insert the UT and diagonal indexes
c
      indexes = 0
      s = 1
      do i = 1, neqns
        o = ptrs(i) + work(i)
             ! o is pointing at the diagonal
        indexes(o) = i
        o = o + 1
        do j = s, s+vss_ptrs(i)-1
             ! j indexes the original structure
          indexes(o) = vss_indexes(j)
          o = o + 1
        end do
        s = j
      end do
c
c           Insert the LT indexes
c
      work = 0
      s = 1
      do c = 1, neqns
        do j = s, s+vss_ptrs(c) - 1
          r = vss_indexes(j)
          ! row r column c, columns will be strictly increasing
          o = ptrs(r) + work(r)
          indexes(o) = c
          work(r) = work(r) + 1
        end do
        s = j
      end do
c
c     All done algorithm is O(nnz), don't think we can do better
c
      return
      end


text:
	PROGRAM TESTLV
C************************************************************************
C* TESTLV								*
C*									*
C* This program tests the LEVEL library subroutines.			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 5/86						*
C* M. desJardins/GSFC	 9/88	Rewritten for GEMPAK4			*
C* M. desJardins/NMC	10/91	Add test for LV_CCRD			*
C* S. Jacobs/NMC	 3/95	Update call to LV_INPT			*
C************************************************************************
	REAL		rlevel (100)
	CHARACTER	vparm*8, clevel*72, vcord*8
	CHARACTER	start*20, stop*20, inc*20
	LOGICAL		mandat
C-------------------------------------------------------------------------
	CALL IN_BDTA ( iret )
	iostat = 0
	DO WHILE  ( iostat .eq. 0 )
	    WRITE  ( 6, 20 )
20	    FORMAT
     +  (/'  1 = LV_CORD   2 = LV_DECD   3 = LV_DUPL   4 = LV_GRNG'/
     +    '  5 = LV_INPT   6 = LV_MANL   7 = LV_SORT   8 = LV_VASL'/
     +    '  9 = LV_CCRD' )
	    CALL TM_INT ( 'Select a subroutine number', .false.,
     +                     .false., 1, numsub, n, ier )
	IF ( ier .eq. 2 ) THEN
           iostat = -1
           numsub = -1
	END IF
C--------------------------------------------------------------------------
	    IF  ( numsub .eq. 1 )  THEN
		WRITE (6,*) 'ENTER VCORD'
		READ  (5,2)  VCORD
		CALL LV_CORD  ( VCORD, VPARM, IVERT, IRET )
		WRITE (6,*) 'VPARM,IVERT,IRET: ',VPARM,IVERT,IRET
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 2 )  THEN
		WRITE (6,*) 'ENTER CLEVEL'
		READ  (5,2)  CLEVEL
2		FORMAT (A)
		CALL LV_DECD  ( CLEVEL, RLEV, IRET )
		WRITE (6,*) 'RLEV,IRET: ',RLEV,IRET
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 3 )  THEN
		WRITE (6,*) 'ENTER NLEV'
		READ  (5,*)  NLEV
		WRITE (6,*) 'ENTER',NLEV,'LEVELS'
		READ  (5,*) (RLEVEL (J),J=1,NLEV)
		CALL LV_DUPL  ( NLEV, RLEVEL, IRET )
		WRITE (6,*) 'NLEV, IRET = ', NLEV, IRET
		WRITE (6,*) 'RLEVEL: ', (RLEVEL(I),I=1,NLEV)
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 4 )  THEN
		WRITE (6,*) 'ENTER START'
		READ  (5,2)  START
		WRITE (6,*) 'ENTER STOP'
		READ  (5,2)  STOP
		WRITE (6,*) 'ENTER INC'
		READ  (5,2)  INC
		WRITE (6,*) 'ENTER NEXP'
		READ  (5,*)  NEXP
		CALL LV_GRNG ( START, STOP, INC, NEXP, RLEVEL, NLEV,
     +				IRET )
		WRITE (6,*) 'NLEV,IRET: ',NLEV, IRET
		WRITE (6,*) 'RLEVEL:',(RLEVEL(I),I=1,NLEV)
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 5 )  THEN
		WRITE (6,*) 'ENTER CLEVEL'
		READ  (5,2)  CLEVEL
		WRITE (6,*) 'ENTER VCOORD'
		READ  (5,2)  VCORD
		WRITE (6,*) 'ENTER NEXP'
		READ  (5,*)  NEXP
		CALL  LV_INPT  ( CLEVEL, NEXP, VCORD, NLEV, RLEVEL,
     +				 LEVTYP, VPARM, IVERT, MANDAT, IRET )
		WRITE (6,*) 'NLEV,LEVTYP,IVERT,MANDAT,IRET',
     +				NLEV,LEVTYP,IVERT,MANDAT,IRET
		IF  ( levtyp .ne. 2 )  THEN
		    WRITE (6,*) 'RLEVEL: ',(RLEVEL(I),I=1,NLEV)
		  ELSE
		    WRITE (6,*) 'RLEVEL: ',(RLEVEL(I),I=1,2)
		END IF
		WRITE (6,*) 'VPARM: ',VPARM
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 6 )  THEN
		WRITE (6,*) 'ENTER NEXP'
		READ  (5,*)  NEXP
		CALL LV_MANL  ( NEXP, NLEV, RLEVEL, IRET )
		WRITE (6,*) 'NLEV,IRET: ',NLEV,IRET
		WRITE (6,*) 'RLEVEL: ',(RLEVEL(I),I=1,NLEV)
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 7 )  THEN
		WRITE (6,*) 'ENTER IVERT,NLEV'
		READ  (5,*)  IVERT,NLEV
		WRITE (6,*) 'ENTER RLEVEL: '
		READ  (5,*)  (RLEVEL(I),I=1,NLEV)
		CALL LV_SORT ( IVERT, NLEV, RLEVEL, IRET )
		WRITE (6,*) 'IRET= ', IRET
		WRITE (6,*) 'RLEVEL: ', (RLEVEL(I),I=1,NLEV)
C------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 8 )  THEN
		WRITE (6,*) 'ENTER NEXP'
		READ  (5,*)  NEXP
		CALL LV_VASL ( NEXP, NLEV, RLEVEL, IRET)
		WRITE (6,*) 'NLEV,IRET: ',NLEV,IRET
		WRITE (6,*) 'RLEVEL: ',(RLEVEL(I),I=1,NLEV)
C-------------------------------------------------------------------------
	      ELSE IF  ( numsub .eq. 9 )  THEN
		WRITE (6,*) 'ENTER IVCORD'
		READ  (5,*)  IVCORD
		CALL LV_CCRD ( IVCORD, VCORD, IRET )
		WRITE (6,*) 'IVCORD, VCOORD, IRET: ', IVCORD, ' ', VCORD,
     +				IRET
C-------------------------------------------------------------------------
	    END IF
	END DO
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
   subroutine foo(func,p,eval)
      real(kind=kind(1.0d0)), dimension(3,0:4,0:4,0:4) :: p
      logical(kind=kind(.true.)), dimension(5,5,5) :: eval
      interface
         subroutine func(values,pt)
            real(kind=kind(1.0d0)), dimension(:), intent(out) :: values
            real(kind=kind(1.0d0)), dimension(:,:), intent(in) :: pt
         end subroutine
      end interface
      real(kind=kind(1.0d0)), dimension(125,3) :: pt
      integer(kind=kind(1)) :: n_pt

      n_pt = 1
      pt(1:n_pt,:) = &
         reshape( &
            pack( &
               transpose(reshape(p,(/3,125/))), &
               spread(reshape(eval,(/125/)),dim=2,ncopies=3)), &
            (/n_pt,3/))

   end subroutine
   end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Describe life here.

Where to begin??

is a series of dogs. George Carlinoon

I was walking my dog at Christmas time in the cloudy afternoon of Christmas Eve in the Village Homes Shire, when I caught a whiff of wood smoke from a fireplace.  That was life.

When I was a little boy, my mother used to sing me a song. It went like this: Life is short, life is shit, and soon it will be over. Cabbie from Brain Candy

Life is like a box of chocolates. A cheap, thoughtless, perfunctory gift that nobody ever asks for. Unreturnable  because all you get back is another box of chocolates. Youre stuck with this undefinable whippedmint crap that you mindlessly wolf down when theres nothing else left to eat. Sure, once in a while, theres a peanut butter cup. An English toffee. But theyre gone too fast, the taste is fleeting. So you end up with nothing but broken bits, filled with hardened jelly and teethshattering nuts. If youre desperate enough to eat those, all youve got left... is an empty box... filled with useless, brown paper wrappers.  The XFiles

Mikey likes it.

  While UCD students can formally study life at the School of Medicine or a biology or philosophy course, you might also consider navigating davis getting out into the world, and adventures outside of davis exploring it, maybe love meeting some people and experiencing life first handcows .

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search