text:
/**
 * @copyright (c) 2017 King Abdullah University of Science and Technology (KAUST).
 *                     All rights reserved.
 **/
/**
 * @file codelet_zgytlr.c
 *
 *  HiCMA codelets kernel
 *  HiCMA is a software package provided by King Abdullah University of Science and Technology (KAUST)
 *
 * @version 0.1.0
 * @author Kadir Akbudak
 * @date 2017-11-16
 * @precisions normal z -> c d s
 **/
#include "morse.h"
#include "runtime/starpu/chameleon_starpu.h"
//#include "runtime/starpu/include/runtime_codelet_z.h"

#include "runtime/starpu/runtime_codelets.h"
ZCODELETS_HEADER(gytlr)

#include <assert.h>
#include <stdio.h>
#include <sys/time.h>//FIXME for gettimeofday
#include <stdlib.h>//FIXME for malloc

#include "starsh.h"
#include "starsh-spatial.h"
#include "starsh-randtlr.h"
#ifdef MKL
  #include <mkl.h>
  #include <mkl_lapack.h>
  //#pragma message("MKL is used")
#else
  #ifdef ARMPL
    #include <armpl.h>
  #else
    #include <cblas.h>
  #endif
  #ifdef LAPACKE_UTILS
    #include <lapacke_utils.h>
  #endif
  #include <lapacke.h>
  //#pragma message("MKL is NOT used")
#endif
extern STARSH_blrf *mpiF;
extern int print_index;
int print_index_end;
extern int store_only_diagonal_tiles;
extern int global_check;
extern int print_mat;
extern void _printmat(double * A, int m, int n, int ld);

//extern int global_always_fixed_rank; //FIXME this does not work, I do not know why
//extern int global_fixed_rank;
int global_always_fixed_rank;
int global_fixed_rank;
int gytlr_tile_ii = -1;
int gytlr_tile_jj = -1;

void zgytlr( int m, int n, /*dimension of squareAD*/
        double *AU,
        double *AV,
        double *AD,
        double *Ark,
        int lda,
        int ldu,
        int ldv,
        int bigM, int ii, int jj, unsigned long long int seed,
        int maxrank, double tol, int compress_diag,
        double *Dense
        )
{
    if(gytlr_tile_ii >= 0) {
        ii = gytlr_tile_ii;
        printf("%s %d: Using fixed i:%d\n", __FILE__, __LINE__, ii);
    }
    if(gytlr_tile_jj >= 0) {
        jj = gytlr_tile_jj;
        printf("%s %d: Using fixed j:%d\n", __FILE__, __LINE__, jj);
    }
    int64_t i, j;
    //printf("m:%d n:%d bigM:%d m0:%d n0:%d\n", m, n, bigM, m0, n0);
    struct timeval tvalBefore, tvalAfter;  // removed comma
    gettimeofday (&tvalBefore, NULL);
    if(print_index){
        fprintf(stderr, "+GYTLR\t|(%d,%d) m:%d n:%d lda:%d ldu:%d ldv:%d\n", ii, jj, m, n, lda, ldu, ldv);
    }

    int shape[2];
    int rank = 0;
    int oversample = 10;
    double *work;
    int *iwork;
    STARSH_cluster *RC = mpiF->row_cluster, *CC = RC;
    void *RD = RC->data, *CD = RD;
    double *saveAD;
                                                    // allocate space for dense tile
    if((ii != jj && store_only_diagonal_tiles == 1) // if tile is off diagonal and
                                                    // and only diagonal tiles are stored in a tall and skinny matrix
                                                    // store_only_diagonal_tiles is here because
                                                    // code can use AD to store dense tiles
            ||                                      // OR
            compress_diag == 1) {                   // diagonals are also compressed so AD may not used perhaps
        saveAD = AD;
        //AD = malloc(sizeof(double) * m * n);
        AD = malloc(sizeof(double) * lda * n);
        assert(m==lda);
    }
    //starsh_blrf_get_block(mpiF, ii, jj, shape, &AD);

    mpiF->problem->kernel(m, n, RC->pivot+RC->start[ii], CC->pivot+CC->start[jj],
            RD, CD, AD, lda);

/*    {*/
        /*if (ii != jj || compress_diag == 1) { */
            /*if(store_only_diagonal_tiles == 1) {*/
                /*assert(AD != saveAD);*/
                /*free(AD);*/
            /*}*/
        /*}*/
        /*return; //TODO*/
    /*}*/
    if(global_check == 1){
       char chall = 'A';
       dlacpy_(&chall, &m, &n, AD, &lda, Dense, &lda);
       //printf("Original problem is copied :%d,%d\n", ii,jj);
    }
    int mn = m;
    int mn2 = maxrank+oversample;
    if(mn2 > mn)
        mn2 = mn;
    // Get size of temporary arrays
    size_t lwork = n, lwork_sdd = (4*mn2+7)*mn2;
    if(lwork_sdd > lwork)
        lwork = lwork_sdd;
    lwork += (size_t)mn2*(2*n+m+mn2+1);
    size_t liwork = 8*mn2;
    // Allocate temporary arrays
    //STARSH_MALLOC(iwork, liwork);
    iwork = malloc(sizeof(*iwork) * liwork);
    if(iwork == NULL) {
        fprintf(stderr, "%s %s %d:\t Allocation failed. No memory! liwork:%d", __FILE__, __func__, __LINE__, liwork);
        exit(-1);
    }
    //STARSH_MALLOC(work, lwork);
    work = malloc(sizeof(*work) * lwork);
    if(work == NULL) {
        fprintf(stderr, "%s %s %d:\t Allocation failed. No memory! lwork:%d", __FILE__, __func__, __LINE__, lwork);
        exit(-1);
    }
    if (ii != jj || compress_diag == 1) { // do not try to compress diagonal blocks if it is not enforced
        //AD is m x n. AU and AV are m x maxrank and n x maxrank correspondingly
        //starsh_kernel_drsdd(m, n, AD, AU, AV, &rank, maxrank, oversample, tol, work, lwork, iwork);
        //starsh_dense_dlrrsdd(m, n, AD, AU, AV, &rank, maxrank, oversample, tol, work, lwork, iwork);
        starsh_dense_dlrrsdd(m, n, AD, lda, AU, ldu,  AV, ldv, &rank, maxrank, oversample, tol, work, lwork, iwork);

    if(0)cblas_dgemm( //for testing purposes
            CblasColMajor,
            CblasNoTrans, CblasTrans,
            m, n, rank,
            1.0,  AU,  ldu,
            AV,  ldv,
            0.0, Dense, lda);

        if(rank == -1){ //means that tile is dense.
            rank = m;
            fprintf(stderr, "%s %s %d: Dense off-diagonal block (%d,%d). maxrank:%d\n", __FILE__, __func__, __LINE__, ii, jj, maxrank);
            exit(0);
        }
        if(rank == 0) rank = 1;
        Ark[0] = rank;
        if(store_only_diagonal_tiles == 1) {
            assert(AD != saveAD);
            free(AD);
        }
        if(global_always_fixed_rank == 1) {
            global_fixed_rank = rank;
        }
        if(print_mat){
            printf("%d\tgytlr-UV-output\n", __LINE__);
            _printmat(AD, m, n, lda);
            _printmat(AU, m, rank, ldu);
            _printmat(AV, ldv, rank, ldv);
        }
    } else {
        Ark[0] = m;
        if(print_mat){
            printf("%d\tgytlr-DENSE-output\n", __LINE__);
            _printmat(AD, m, m, lda);
        }
    }

    /*printf("m:%d n:%d Tile %d,%d rk:%d lda:%d maxrank:%d oversample:%d tol:%.2e AD:%p AU:%p AV:%p\n", */
            /*m, n, ii, jj, rank, lda, maxrank, oversample, tol, AD, AU, AV);*/
    /*double *tmp = AD;*/
    /*for (i = 0; i < m; ++i) {*/
        /*for (j=0; j<n; ++j ) {*/
            /*printf("%+.2e ",tmp[j*lda+i]);*/
        /*}*/
        /*printf("\n");*/
    /*}*/
    /*printf("\n");*/

    free(work);
    free(iwork);
    if(print_index || print_index_end ){
        gettimeofday (&tvalAfter, NULL);
        fprintf(stderr, "-GYTLR\t|(%d,%d) rk:%g m:%d n:%d lda:%d ldu:%d ldv:%d\t\t\t\t\tGYTLR: %.4f\n",ii,jj,Ark[0],m, n, lda, ldu, ldv,
                (tvalAfter.tv_sec - tvalBefore.tv_sec)
                 +(tvalAfter.tv_usec - tvalBefore.tv_usec)/1000000.0
                );
    }
}

/*   MORSE_TASK_zgytlr - Generate a tile for random matrix. */

void HICMA_TASK_zgytlr( const MORSE_option_t *options,
                        int m, int n,
                        const MORSE_desc_t *AUV,
                        const MORSE_desc_t *Ark,
                        int Am, int An,
                        int lda,
                        int ldu,
                        int ldv,
                        int bigM, int m0, int n0, unsigned long long int seed,
                        int maxrank, double tol,
                        int compress_diag,
                        MORSE_desc_t *Dense
                        )
{
    struct starpu_codelet *codelet = &cl_zgytlr;
    //void (*callback)(void*) = options->profiling ? cl_zgytlr_callback : NULL;
    void (*callback)(void*) = NULL;
    int nAUV = AUV->nb;

    MORSE_BEGIN_ACCESS_DECLARATION;
    MORSE_ACCESS_W(AUV, Am, An);
    MORSE_ACCESS_W(Ark, Am, An);
    MORSE_ACCESS_RW(Dense, Am, An);
    MORSE_END_ACCESS_DECLARATION;

    //printf("%s:%d: Am:%d An:%d lda:%d bigM:%d m0:%d n0:%d\n ", __FILE__, __LINE__, Am, An, lda, bigM, m0, n0);
        //printf("%s %d: Am:%d An:%d ADm:%d ADn:%d ptr:%p\n", __func__, __LINE__, Am, An, ADm, ADn, ptr);
    starpu_insert_task(
            starpu_mpi_codelet(codelet),
            STARPU_VALUE,    &m,                      sizeof(int),
            STARPU_VALUE,    &n,                      sizeof(int),
            STARPU_VALUE,    &nAUV,                      sizeof(int),
            STARPU_W,         RTBLKADDR(AUV, double, Am, An),
            STARPU_W,         RTBLKADDR(Ark, double, Am, An),
            STARPU_RW,         RTBLKADDR(Dense, double, Am, An), // _R must be _W SERIOUSLY. BUT _W STALLS ON SHAHEEN. FIXME
            STARPU_VALUE,  &lda,                      sizeof(int),
            STARPU_VALUE,  &ldu,                      sizeof(int),
            STARPU_VALUE,  &ldv,                      sizeof(int),
            STARPU_VALUE, &bigM,                      sizeof(int),
            STARPU_VALUE,   &Am,                      sizeof(int),
            STARPU_VALUE,   &An,                      sizeof(int),
            STARPU_VALUE, &seed,   sizeof(unsigned long long int),
            STARPU_VALUE,   &maxrank,                      sizeof(int),
            STARPU_VALUE,   &tol,                      sizeof(double),
            STARPU_VALUE,   &compress_diag,            sizeof(int),
            STARPU_PRIORITY,    options->priority,
            STARPU_CALLBACK,    callback,
#if defined(CHAMELEON_CODELETS_HAVE_NAME)
            STARPU_NAME, "zgytlr",
#endif
            0);
}

/*   cl_zgytlr_cpu_func - Generate a tile for random matrix. */

#if !defined(CHAMELEON_SIMULATION)
static void cl_zgytlr_cpu_func(void *descr[], void *cl_arg)
{
    int m;
    int n;
    int nAUV;
    double *AUV;
    double *AD = NULL;
    double *Ark;
    double *Dense;
    int lda;
    int ldu;
    int ldv;
    int bigM;
    int m0;
    int n0;
    unsigned long long int seed;
    int maxrank;
    double tol;
    int compress_diag;

    AUV = (double *)STARPU_MATRIX_GET_PTR(descr[0]);
    Ark = (double *)STARPU_MATRIX_GET_PTR(descr[1]);
    Dense = (double *)STARPU_MATRIX_GET_PTR(descr[2]);

    starpu_codelet_unpack_args(cl_arg, &m, &n, &nAUV, &lda, &ldu, &ldv, &bigM, &m0, &n0, &seed, &maxrank, &tol, &compress_diag );

    double *AU = AUV;
    int nAU = nAUV/2;
    assert(ldu == ldv);
    size_t nelm_AU = (size_t)ldu * (size_t)nAU;
    double *AV = &(AUV[nelm_AU]);

    //printf("(%d,%d)%d %s %d %d\n", m0/m,n0/n,MORSE_My_Mpi_Rank(), __func__, __LINE__, AD == Dense);
    zgytlr( m, n,
            AU,
            AV,
            AD,
            Ark,
            lda,
            ldu,
            ldv,
            bigM, m0, n0, seed,
            maxrank, tol,
            compress_diag,
            Dense
            );
}
#endif /* !defined(CHAMELEON_SIMULATION) */

/*
 * Codelet definition
 */
CODELETS_CPU(zgytlr, 3, cl_zgytlr_cpu_func)


text:
#ifndef IAIF_H
#define IAIF_H

#include <utility>
#include <tuple>

#include <gsl/gsl_vector.h>

#include "constants.h"
#include "lpc.h"
#include "filter.h"

void computeIAIF(gsl_vector *g, gsl_vector *dg, gsl_vector *x);

#endif // IAIF_H

ext: h

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

BC_INLINE2(GSL_MONTE_FN_EVAL,gsl_monte_function*,double*,double)

ext: c


text:
/**
 *
 * @file core_dssssm.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:45 2014
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dssssm applies the LU factorization update from a complex
 *  matrix formed by a lower triangular IB-by-K tile L1 on top of a
 *  M2-by-K tile L2 to a second complex matrix formed by a M1-by-N1
 *  tile A1 on top of a M2-by-N2 tile A2 (N1 == N2).
 *
 *  This is the right-looking Level 2.5 BLAS version of the algorithm.
 *
 *******************************************************************************
 *
 * @param[in] M1
 *         The number of rows of the tile A1.  M1 >= 0.
 *
 * @param[in] N1
 *         The number of columns of the tile A1.  N1 >= 0.
 *
 * @param[in] M2
 *         The number of rows of the tile A2 and of the tile L2.
 *         M2 >= 0.
 *
 * @param[in] N2
 *         The number of columns of the tile A2.  N2 >= 0.
 *
 * @param[in] K
 *         The number of columns of the tiles L1 and L2.  K >= 0.
 *
 * @param[in] IB
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in,out] A1
 *         On entry, the M1-by-N1 tile A1.
 *         On exit, A1 is updated by the application of L (L1 L2).
 *
 * @param[in] LDA1
 *         The leading dimension of the array A1.  LDA1 >= max(1,M1).
 *
 * @param[in,out] A2
 *         On entry, the M2-by-N2 tile A2.
 *         On exit, A2 is updated by the application of L (L1 L2).
 *
 * @param[in] LDA2
 *         The leading dimension of the array A2.  LDA2 >= max(1,M2).
 *
 * @param[in] L1
 *         The IB-by-K lower triangular tile as returned by
 *         CORE_dtstrf.
 *
 * @param[in] LDL1
 *         The leading dimension of the array L1.  LDL1 >= max(1,IB).
 *
 * @param[in] L2
 *         The M2-by-K tile as returned by CORE_dtstrf.
 *
 * @param[in] LDL2
 *         The leading dimension of the array L2.  LDL2 >= max(1,M2).
 *
 * @param[in] IPIV
 *         The pivot indices array of size K as returned by
 *         CORE_dtstrf.
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dssssm = PCORE_dssssm
#define CORE_dssssm PCORE_dssssm
#endif
int CORE_dssssm(int M1, int N1, int M2, int N2, int K, int IB,
                double *A1, int LDA1,
                double *A2, int LDA2,
                const double *L1, int LDL1,
                const double *L2, int LDL2,
                const int *IPIV)
{
    static double zone  = 1.0;
    static double mzone =-1.0;

    int i, ii, sb;
    int im, ip;

    /* Check input arguments */
    if (M1 < 0) {
        coreblas_error(1, "Illegal value of M1");
        return -1;
    }
    if (N1 < 0) {
        coreblas_error(2, "Illegal value of N1");
        return -2;
    }
    if (M2 < 0) {
        coreblas_error(3, "Illegal value of M2");
        return -3;
    }
    if (N2 < 0) {
        coreblas_error(4, "Illegal value of N2");
        return -4;
    }
    if (K < 0) {
        coreblas_error(5, "Illegal value of K");
        return -5;
    }
    if (IB < 0) {
        coreblas_error(6, "Illegal value of IB");
        return -6;
    }
    if (LDA1 < max(1,M1)) {
        coreblas_error(8, "Illegal value of LDA1");
        return -8;
    }
    if (LDA2 < max(1,M2)) {
        coreblas_error(10, "Illegal value of LDA2");
        return -10;
    }
    if (LDL1 < max(1,IB)) {
        coreblas_error(12, "Illegal value of LDL1");
        return -12;
    }
    if (LDL2 < max(1,M2)) {
        coreblas_error(14, "Illegal value of LDL2");
        return -14;
    }

    /* Quick return */
    if ((M1 == 0) || (N1 == 0) || (M2 == 0) || (N2 == 0) || (K == 0) || (IB == 0))
        return PLASMA_SUCCESS;

    ip = 0;

    for(ii = 0; ii < K; ii += IB) {
        sb = min(K-ii, IB);

        for(i = 0; i < sb; i++) {
            im = IPIV[ip]-1;

            if (im != (ii+i)) {
                im = im - M1;
                cblas_dswap(N1, &A1[ii+i], LDA1, &A2[im], LDA2);
            }
            ip = ip + 1;
        }

        cblas_dtrsm(
            CblasColMajor, CblasLeft, CblasLower,
            CblasNoTrans, CblasUnit,
            sb, N1, (zone),
            &L1[LDL1*ii], LDL1,
            &A1[ii], LDA1);

        cblas_dgemm(
            CblasColMajor, CblasNoTrans, CblasNoTrans,
            M2, N2, sb,
            (mzone), &L2[LDL2*ii], LDL2,
            &A1[ii], LDA1,
            (zone), A2, LDA2);
    }
    return PLASMA_SUCCESS;
}


text:
/*
  Copyright [2017-2020] [IBM Corporation]
  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
  http://www.apache.org/licenses/LICENSE-2.0
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
*/
#ifndef __SHARD_CONNECTION_H__
#define __SHARD_CONNECTION_H__

#ifdef __cplusplus

#include "tls_session.h"
#include "buffer_manager.h"
#include "fabric_connection_base.h"  // default to fabric transport
#include "mcas_config.h"
#include "pool_manager.h"
#include "protocol.h"
#include "protocol_ostream.h"
#include "region_manager.h"
#include "connection_state.h"

#include <api/components.h>
#include <api/fabric_itf.h>
#include <api/kvstore_itf.h>
#include <common/exceptions.h>
#include <common/logging.h>
#include <gsl/pointers>

#include <cassert>
#include <map>
#include <queue>
#include <utility> /* swap */
#include <vector>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"

namespace mcas
{
using Connection_base = Fabric_connection_base;

/**
 * Connection handler is instantiated for each "connected" client
 */
class Connection_handler
  : public Connection_base,
    public Region_manager,
    private Connection_TLS_session
{
  friend class Connection_TLS_session;
  friend struct TLS_transport;

public:
  enum {
        TICK_RESPONSE_CONTINUE        = 0,
        TICK_RESPONSE_BOOTSTRAP_SPAWN = 1,
        TICK_RESPONSE_WAIT_SECURITY   = 2,
        TICK_RESPONSE_CLOSE           = 0xFF,
  };

  enum {
        ACTION_NONE = 0,
        ACTION_RELEASE_VALUE_LOCK_EXCLUSIVE,
        ACTION_POOL_DELETE,
  };

private:
  Connection_state    _state       = Connection_state::INITIAL;
  unsigned option_DEBUG = mcas::global::debug_level;

  /* list of pre-registered memory regions; normally one region */
  std::vector<component::IKVStore::memory_handle_t> _mr_vector;

  uint64_t                            _tick_count alignas(8);
  uint64_t                            _auth_id;
  std::queue<buffer_t *>              _pending_msgs;
  std::queue<action_t>                _pending_actions;
  Pool_manager                        _pool_manager; /* instance shared across connections */

  /* Adaptor point for different transports */
  using Connection = component::IFabric_server;
  using Factory    = component::IFabric_server_factory;

  struct stats {
    uint64_t response_count;
    uint64_t recv_msg_count;
    uint64_t send_msg_count;
    uint64_t wait_send_value_misses;
    uint64_t wait_msg_recv_misses;
    uint64_t wait_respond_complete_misses;
    uint64_t last_count;
    uint64_t next_stamp;
    stats()
      : response_count(0),
        recv_msg_count(0),
        send_msg_count(0),
        wait_send_value_misses(0),
        wait_msg_recv_misses(0),
        wait_respond_complete_misses(0),
        last_count(0),
        next_stamp(0)
    {
    }
  } _stats alignas(8);

  void dump_stats()
  {
    PINF("-----------------------------------------");
    PINF("| Connection Handler Statistics         |");
    PINF("-----------------------------------------");
    PINF("Ticks                       : %lu", _tick_count);
    PINF("Open pools                  : %lu", _pool_manager.open_pool_count());
    PINF("NEW_MSG_RECV misses         : %lu", _stats.wait_msg_recv_misses);
    PINF("Recv message count          : %lu", _stats.recv_msg_count);
    PINF("Send message count          : %lu", _stats.send_msg_count);
    PINF("Response count              : %lu", _stats.response_count);
    PINF("WAIT_SEND_VALUE misses      : %lu", _stats.wait_send_value_misses);
    PINF("WAIT_RESPOND_COMPLETE misses: %lu", _stats.wait_respond_complete_misses);
    PINF("-----------------------------------------");
  }

  /**
   * Change state in FSM
   *
   * @param s
   */
  inline void set_state(Connection_state s)
  {
    if (2 < option_DEBUG) {
      const std::map<Connection_state, const char *> m{
                                          {Connection_state::INITIAL, "INITIAL"},
                                          {Connection_state::WAIT_HANDSHAKE, "WAIT_HANDSHAKE"},
                                          {Connection_state::WAIT_TLS_HANDSHAKE, "WAIT_TLS_HANDSHAKE"},
                                          {Connection_state::CLIENT_DISCONNECTED, "CLIENT_DISCONNECTED"},
                                          {Connection_state::CLOSE_CONNECTION, "CLOSE_CONNECTION"},
                                          {Connection_state::WAIT_NEW_MSG_RECV, "WAIT_NEW_MSG_RECV"}};
      PLOG("state %s -> %s", m.find(_state)->second, m.find(s)->second);
    }
    _state = s; /* we could add transition checking later */
  }

  static void static_send_callback2(void *cnxn, buffer_t *iob) noexcept
  {
    auto base = static_cast<Fabric_connection_base *>(cnxn);
    static_cast<Connection_handler *>(base)->send_callback2(iob);
  }

  void send_callback2(buffer_t *iob) noexcept
  {
    assert(iob->value_adjunct);
    if (0 < option_DEBUG) {
      PLOG("Completed send2 (value_adjunct %p)", common::p_fmt(iob->value_adjunct));
    }
    _deferred_unlock.push(iob->value_adjunct);
    send_callback(iob);
  }

  /**
   * Send handshake response
   *
   * @param start_tls Set if a TLS handshake needs to be started
   */
  void respond_to_handshake(bool start_tls);

public:

  explicit Connection_handler(unsigned debug_level,
                              gsl::not_null<Factory *> factory,
                              gsl::not_null<Connection *> connection);

  virtual ~Connection_handler();

  inline bool client_connected() { return _state != Connection_state::CLIENT_DISCONNECTED; }

  auto allocate_send() { return allocate(static_send_callback); }
  auto allocate_recv() { return allocate(static_recv_callback); }

  /**
   * Initialize security session
   *
   * @param bind_ipaddr
   * @param bind_port
   */
  void configure_security(const std::string& bind_ipaddr,
                          const unsigned bind_port,
                          const std::string& cert_file,
                          const std::string& key_file);

  /**
   * State machine transition tick.  It is really important that this tick
   * execution duration is small, so that other connections are not impacted by
   * the thread not returning to them.
   *
   */
  int tick();
  /**
   * Check for network completions
   *
   */
  Fabric_connection_base::Completion_state check_network_completions()
  {
    auto state = poll_completions();

    while (!_deferred_unlock.empty()) {
      //      void *deferred_unlock = nullptr;
      //???     std::swap(deferred_unlock, _deferred_unlock.front());
      void *deferred_unlock = _deferred_unlock.front();
      if (option_DEBUG > 2)
        PLOG("adding action for deferred unlocking value @ %p", deferred_unlock);
      add_pending_action(action_t{ACTION_RELEASE_VALUE_LOCK_EXCLUSIVE, deferred_unlock});
      _deferred_unlock.pop();
    }

    return state;
  }

  /**
   * Peek at pending message from the connection. Used when the resources
   * required to process a pending nessage may depend on the content of the
   * message.
   *
   * @param msg [out] Pointer to base protocol message
   *
   * @return Pointer to buffer holding the message or null if there are none
   */
  inline mcas::protocol::Message *peek_pending_msg() const
  {
    return _pending_msgs.empty() ? nullptr
      : static_cast<mcas::protocol::Message *>(_pending_msgs.front()->base().get());
  }

  /**
   * Discard a pending message from the connection. Used along with tick
   * (which adds pending messages) and peek (which peeks at them).
   *
   * @param msg [out] Pointer to base protocol message
   *
   * @return Pointer to buffer holding the message or null if there are none
   */
  inline buffer_t *pop_pending_msg()
  {
    assert(!_pending_msgs.empty());
    auto iob = _pending_msgs.front();
    _pending_msgs.pop();
    return iob;
  }

  /**
   * Get deferred actions
   *
   * @param action [out] Action
   *
   * @return True if action
   */
  inline bool get_pending_action(action_t &action)
  {
    if (_pending_actions.empty()) return false;

    action = _pending_actions.front();

    if (option_DEBUG > 2) PLOG("Connection_handler: popped pending action (%u, %p)", action.op, action.parm);

    _pending_actions.pop();
    return true;
  }

  template <typename MT>
  void msg_log(const unsigned level_, const MT *msg_, const char *desc_, const char *direction_)
  {
    if (level_ < Connection_base::debug_level()) {
      std::ostringstream m;
      m << *msg_;
      PLOG("%s (%s) %s", direction_, desc_, m.str().c_str());
    }
  }

  template <typename MT>
  void msg_recv_log(const MT *m, const char *desc)
  {
    msg_recv_log(1, m, desc);
  }

  template <typename MT>
  void msg_recv_log(const unsigned level, const MT *msg, const char *desc)
  {
    msg_log(level, msg, desc, "RECV");
  }

  template <typename MT>
  void msg_send_log(const MT *m, const char *desc)
  {
    msg_send_log(1, m, desc);
  }

  template <typename MT>
  void msg_send_log(const unsigned level, const MT *msg, const char *desc)
  {
    msg_log(level, msg, desc, "SEND");
  }

  /**
   * Add an action to the pending queue
   *
   * @param action Action to add
   */
  inline void add_pending_action(const action_t& action)
  {
    _pending_actions.push(action);
  }

  template <typename Msg>
  void post_send_buffer(gsl::not_null<buffer_t *> buffer, Msg *msg, const char *desc)
  {
    msg_send_log(msg, desc);
    Connection_base::post_send_buffer(buffer);
  }

  template <typename Msg>
  void post_send_buffer_with_value(gsl::not_null<buffer_t *> buffer,
                         const ::iovec &           val_iov,
                         void *                    val_desc,
                         Msg *                     msg,
                         const char *              func_name)
  {
    msg_send_log(msg, func_name);
    Connection_base::post_send_buffer2(buffer, val_iov, val_desc);
  }

  /**
   * Post a response
   *
   * @param iob IO buffer to post
   */
  template <typename Msg>
  inline void post_response_with_value(gsl::not_null<buffer_t *> iob,
                             const ::iovec &           val_iov,
                             void *                    val_desc,
                             Msg *                     msg,
                             const char *              func_name)
  {
    iob->set_completion(static_send_callback2, val_iov.iov_base);
    post_send_buffer_with_value(iob, val_iov, val_desc, msg, func_name);
    /* don't wait for this, let it be picked up in the check_completions cycle
     */

    _stats.response_count++;
  }

  /**
   * Post a response
   *
   * @param iob IO buffer to post
   */
  template <typename Msg>
  inline void post_response(gsl::not_null<buffer_t *> iob, Msg *msg, const char *desc)
  {
    post_send_buffer(iob, msg, desc);
    /* don't wait for this, let it be picked up in the check_completions cycle
     */

    _stats.response_count++;
  }

  /**
   * Set up for pending value send/recv
   *
   * @param target
   * @param target_len
   * @param region
   */
  void add_memory_handle(component::IKVStore::memory_handle_t handle) { _mr_vector.push_back(handle); }

  component::IKVStore::memory_handle_t pop_memory_handle()
  {
    if (_mr_vector.empty()) return nullptr;
    auto mr = _mr_vector.back();
    _mr_vector.pop_back();
    return mr;
  }

  inline uint64_t       auth_id() const { return _auth_id; }
  inline void           set_auth_id(uint64_t id) { _auth_id = id; }
  inline size_t         max_message_size() const { return _max_message_size; }
  inline Pool_manager & pool_manager() { return _pool_manager; }

private:
  common::Byte_buffer _tls_buffer;
};

}  // namespace mcas

#pragma GCC diagnostic pop

#endif

#endif  // __CONNECTION_HANDLER_H__


text:
#include <stdio.h>
#include <stdlib.h>
#include "parmt_utils.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "iscl/array/array.h"

/*!
 * @brief The iloc'th location and iobs'th observation this sets the forward
 *        modeling matrix G s.t. Gm = u
 *
 * @param[in] iobs    observation number
 * @param[in] iloc    location number
 * @param[in] npmax   leading dimension of G
 * @param[in] data    data data structure with Green's functions which will be
 *                    set on G
 *
 * @param[out] G      column major modeling matrix s.t. Gm=u where
 *                    G is packed
 *                     \f$ \{ G_{xx}, G_{yy}, G_{zz},
 *                            G_{xy}, G_{xz}, G_{yz} \} \f$
 *                    and m is in NED ordering [npmax x 6]
 *
 * @result 0 indicates success
 *
 */
int parmt_utils_setDataOnG(const int iobs, const int iloc, const int npmax,
                           const struct parmtData_struct data,
                           double *__restrict__ G)
{
    const char *fcnm = "parmt_utils_setDataOnG\0";
    int k;
    if (iobs < 0 || iobs >= data.nobs)
    {
        printf("%s: Error invalid observation\n", fcnm);
        return -1;
    }
    if (iloc < 0 || iloc >= data.nlocs)
    {
        printf("%s: Error invalid location\n", fcnm);
        return -1;
    }
    k = iobs*data.nlocs + iloc;
    if (npmax > data.sacGxx[k].npts || npmax > data.sacGyy[k].npts ||
        npmax > data.sacGzz[k].npts || npmax > data.sacGxy[k].npts ||
        npmax > data.sacGxz[k].npts || npmax > data.sacGyz[k].npts)
    {
        printf("%s: npmax is too small\n", fcnm);
        return -1;
    }
    if (data.sacGxx[k].data == NULL || data.sacGyy[k].data == NULL ||
        data.sacGzz[k].data == NULL || data.sacGxy[k].data == NULL ||
        data.sacGxz[k].data == NULL || data.sacGyz[k].data == NULL)
    {
        printf("%s: Greens function is null\n", fcnm);
        return -1;
    }
    // Get the Green's functions onto the matrix G
    array_zeros64f_work(6*npmax, G);
    cblas_dcopy(data.sacGxx[k].header.npts,
                data.sacGxx[k].data, 1, &G[0*npmax], 1);
    cblas_dcopy(data.sacGyy[k].header.npts,
                data.sacGyy[k].data, 1, &G[1*npmax], 1);
    cblas_dcopy(data.sacGzz[k].header.npts,
                data.sacGzz[k].data, 1, &G[2*npmax], 1);
    cblas_dcopy(data.sacGxy[k].header.npts,
                data.sacGxy[k].data, 1, &G[3*npmax], 1);
    cblas_dcopy(data.sacGxz[k].header.npts,
                data.sacGxz[k].data, 1, &G[4*npmax], 1);
    cblas_dcopy(data.sacGyz[k].header.npts,
                data.sacGyz[k].data, 1, &G[5*npmax], 1);
    return 0;
}


text:
/*
 * Tema 2 ASC
 * 2019 Spring
 * Catalin Olaru / Vlad Spoiala
 */
#include <cblas.h>

#include "utils.h"
/*
 * Add your BLAS implementation here
 */
double* my_solver(int N, double *A, double *B) {
	double *first_multiplication = malloc(N * N * sizeof(double));
	double *second_multiplication = malloc(N * N * sizeof(double));
	double *final_result = malloc(N * N * sizeof(double));
	double alpha = 1.0;
	double beta = 0.0;
	double *identity_matrix = malloc(N * N * sizeof(double));
	int i, j;

/* Compute identity matrix so we can use it later on. */
	for (i = 0; i < N; ++i)
		for (j = 0; j < N; ++j)
			identity_matrix[i * N + j] = j == i ? 1.0 : 0.0;

/* Compute Transpose(A) * B in first_multiplication. */
	cblas_dgemm(CblasRowMajor,
				CblasTrans,
				CblasNoTrans,
				N,
				N,
				N,
				alpha,
				A,
				N,
				B,
				N,
				beta,
				first_multiplication,
				N);
/* Compute Transpose(B) * A in second_multiplication. */
	cblas_dgemm(CblasRowMajor,
				CblasTrans,
				CblasNoTrans,
				N,
				N,
				N,
				alpha,
				B,
				N,
				A,
				N,
				beta,
				second_multiplication,
				N);
/* Add the previous 2 results in second multiplication. */
	cblas_dgemm(CblasRowMajor,
				CblasNoTrans,
				CblasNoTrans,
				N,
				N,
				N,
				alpha,
				first_multiplication,
				N,
				identity_matrix,
				N,
				alpha,
				second_multiplication,
				N);

/* Apply zerotr(second_multiplication) as everything involving
 * CblasLower or CblasUpper didn't work for some reason. :(
 */
	for (i = 0; i < N; ++i)
		for (j = 0; j < i; ++j)
			second_multiplication[i * N + j] = 0;

/* Multiply second_multiplication by itself so we get the
 * result in final_result. */
	cblas_dgemm(CblasRowMajor,
				CblasNoTrans,
				CblasNoTrans,
				N,
				N,
				N,
				alpha,
				second_multiplication,
				N,
				second_multiplication,
				N,
				beta,
				final_result,
				N);

/* Free everything aside from the return value. */
	free(first_multiplication);
	free(second_multiplication);
	free(identity_matrix);

	return final_result;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <math.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_voxel_graphics.h"
#include "tz_stack_code.h"
#include "tz_stack_math.h"
#include "tz_stack_utils.h"
#include "tz_objdetect.h"
#include "tz_voxeltrans.h"
#include "tz_stack_stat.h"
#include "tz_stack_draw.h"
#include "tz_geo3d_vector.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_bwdist.h"
#include "tz_neurotrace.h"
#include "tz_locseg_chain.h"
#include "tz_trace_utils.h"
#include "tz_stack_threshold.h"
#include "tz_darray.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];

#if 1
  sprintf(file_path, "../data/%s/bundle_seed.tif", neuron_name);
  Stack *stack = Read_Stack(file_path);
  Voxel_List *list = Stack_To_Voxel_List(stack);

  sprintf(file_path, "../data/%s/bundle_seeds.pa", neuron_name);
  Pixel_Array *pa = Pixel_Array_Read(file_path);

  Voxel *seed;

  int i;
  double *pa_array = (double *) pa->array;

  printf("mean: %g, std: %g\n", gsl_stats_mean(pa_array, 1, pa->size),
	 sqrt(gsl_stats_variance(pa_array, 1, pa->size)));

  double threshold = gsl_stats_mean(pa_array, 1, pa->size) +
    3.0 * sqrt(gsl_stats_variance(pa_array, 1, pa->size));

  double max_r = gsl_stats_max(pa_array, 1, pa->size);

  printf("%g, %g\n", threshold, max_r);

  max_r = max_r * 2;
  Set_Neuroseg_Max_Radius(max_r / 2.0);

  dim_type dim[3];
  dim[0] = stack->width;
  dim[1] = stack->height;
  dim[2] = stack->depth;

  Kill_Stack(stack);

  sprintf(file_path, "../data/%s/bundle.tif", neuron_name);
  stack = Read_Stack(file_path);

  Rgb_Color color;
  Set_Color(&color, 255, 0, 0);

  sprintf(file_path, "../data/%s.tif", neuron_name);
  Stack *signal = Read_Stack(file_path);

  Stack *canvas = NULL;

  char trace_mask_path[100];
  sprintf(trace_mask_path, "../data/%s/trace_mask.tif", neuron_name);
  char trace_file_path[100];
  sprintf(trace_file_path, "../data/%s/traced.tif", neuron_name);

  if (fexist(trace_file_path) == 1) {
    canvas = Read_Stack((char *) trace_file_path);
  } else {
    canvas = Translate_Stack(signal, COLOR, 0);
  }

  Stack *traced = NULL;

  if (fexist(trace_mask_path) == 1) {
    traced = Read_Stack((char *) trace_mask_path);
    if (traced->kind != GREY16) {
      Translate_Stack(traced, GREY16, 1);
    }
  } else {
    traced = Make_Stack(GREY16, signal->width, signal->height, signal->depth);
    Zero_Stack(traced);
  }

  Neurochain *chain = NULL;

  double z_scale = 1.0;
  sprintf(file_path, "../data/%s.res", neuron_name);
  if (fexist(file_path)) {
    double res[3];
    int length;
    darray_read2(file_path, res, &length);
    if (res[0] != res[1]) {
      perror("Different X-Y resolutions.");
      TZ_ERROR(ERROR_DATA_VALUE);
    }
    z_scale = res[0] / res[2];
  }

  /*
  sprintf(file_path, "../data/%s/soma.tif", neuron_name);
  Stack *soma = Read_Stack(file_path);
  Stack_Threshold(soma, 1);
  Stack_Binarize(soma);
  Stack_Not(soma, soma);
  Stack_And(traced, soma, traced);
  Kill_Stack(soma);
  */

  tic();

  FILE *fp = NULL;
  char chain_file_path[100];
  char vrml_file_path[100];
  int seg_idx = 229;
  int min_chain_length = 4;

#define TRY_NEXT_SEED(current_seed)			\
  Kill_Voxel(current_seed);				\
  continue;

#define STOP_TRACING(current_seed)			\
  Kill_Voxel(current_seed);				\
  break;

  Trace_Workspace *tw = New_Trace_Workspace();
  tw->length = 200;
  tw->fit_first = FALSE;
  tw->tscore_option = 1;
  tw->min_score = 0.4;
  tw->trace_direction = DL_BOTHDIR;
  tw->trace_mask = traced;
  tw->dyvar[0] = max_r;
  tw->test_func = Locseg_Chain_Trace_Test;

  int chain_id = 100;

  char chain_info_file_path[100];
  sprintf(chain_info_file_path, "../data/%s/chain_info.txt", neuron_name);
  FILE *chain_info_file = fopen(chain_info_file_path, "w");

  for (i = 0; i < pa->size; i++) {
    /* do not skip seed dequeueing */
    seed = Voxel_Queue_De(&list);
    printf("---------------------------------> seed: %d / %d\n", i, pa->size);

    sprintf(chain_file_path, "../data/%s/chain%d.bn", neuron_name, i);
    sprintf(vrml_file_path, "../data/%s/chain%d.wrl", neuron_name, i);

    if (fexist(chain_file_path) == 1) {
      chain = Read_Neurochain(chain_file_path);
      if (Neurochain_Length(chain, FORWARD) >= min_chain_length) {
	Write_Neurochain_Vrml(vrml_file_path, chain);
      }
      Neurochain_Label(canvas, chain, z_scale);
      Free_Neurochain(chain);
      TRY_NEXT_SEED(seed);
    }

    /* for debugging */
    if (i < seg_idx) {
      //TRY_NEXT_SEED(seed);
    }
    /****************/

    if (*STACK_PIXEL_16(traced, seed->x, seed->y, seed->z, 0) > 0) {
      TRY_NEXT_SEED(seed);
    }

    double width = pa_array[i];

    if (width > max_r) {
      TRY_NEXT_SEED(seed);
    }

    Print_Voxel(seed);
    printf("%g\n", width);

    int max_level = (int) (width + 0.5);
    if (max_level < 1.1) {
      TRY_NEXT_SEED(seed);
      //max_level = 6;
    }

    /* for debugging */
    if (i > seg_idx) {
      //STOP_TRACING(seed);
    }
    /****************/

    double cpos[3];
    cpos[0] = seed->x;
    cpos[1] = seed->y;
    cpos[2] = seed->z;
    cpos[2] *= z_scale;

    Local_Neuroseg *startseg = New_Local_Neuroseg();;

    Set_Neuroseg_Position(startseg, cpos, NEUROSEG_CENTER);
    Reset_Neuroseg(&(startseg->seg));
    startseg->seg.r1 = width;
    startseg->seg.r2 = width;

    Stack_Fit_Score fs;
    fs.n = 1;
    fs.options[0] = tw->tscore_option;

    printf("Search orientation ...\n");
    Local_Neuroseg_Orientation_Search_C(startseg, signal, z_scale, &fs);

    Print_Local_Neuroseg(startseg);

    printf("Start: \n");

    Locseg_Chain* locseg_chain =
      Locseg_Chain_Trace_Init(signal, z_scale, startseg, &fs);

    if (fs.scores[0] >= tw->min_score) {
      Locseg_Chain_Iterator_Start(locseg_chain, DL_HEAD);
      Local_Neuroseg *head = Locseg_Chain_Next(locseg_chain);
      if ((head->seg.r1 < max_r) && (head->seg.r2 < max_r)) {
	Trace_Locseg(signal, z_scale, locseg_chain, tw);
	chain = Neurochain_From_Locseg_Chain(locseg_chain);
	//chain = Trace_Neuron2(signal, chain, BOTH, traced, z_scale, 100);

	printf("Remove ends ...\n");
	chain = Neurochain_Remove_Overlap_Ends(chain);
	chain = Neurochain_Remove_Turn_Ends(chain, 1.0);

	//Print_Neurochain(Neurochain_Head(chain));
	//Stack_Draw_Object_Bwc(canvas, obj, color);
	Neurochain *chain_head = Neurochain_Head(chain);

	printf("Saving ...\n");
	fp = fopen(chain_file_path, "w");
	Neurochain_Fwrite(chain_head, fp);
	fclose(fp);
	if (Neurochain_Length(chain_head, FORWARD) >= min_chain_length) {
	  printf("Labeling ...\n");
	  Write_Neurochain_Vrml(vrml_file_path, chain_head);
	  Neurochain_Label_G(trac, chain_head, z_scale, 0,
			     Neurochain_Length(chain_head, FORWARD),
			     1.5, 0.0, chain_id);
	  fprintf(chain_info_file, "%d %d %d\n", i, chain_id,
		  Neurochain_Length(chain_head, FORWARD));
	  chain_id++;
	  Neurochain_Label(canvas, chain_head, z_scale);
	  Write_Stack((char *) trace_mask_path, traced);
	}
	//Neurochain_Erase(traced, Neurochain_Head(chain), z_scale);
	Free_Neurochain(chain);
      }
    } else {
      printf("Ignored seed:\n");

    }

    Kill_Locseg_Chain(locseg_chain);
    Kill_Voxel(seed);
  }

  fclose(chain_info_file);

  Write_Stack((char *) trace_file_path, canvas);

  Kill_Pixel_Array(pa);
  Voxel_List_Removeall(&list);

  printf("Time passed: %lld\n", toc());

#endif

  return 0;
}


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <mpi.h>
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
#include <petscpc.h>
#include <petscsnes.h>

#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include <petsc/private/kspimpl.h>
  #else
     #include <petsc-private/kspimpl.h>   /*I "petscksp.h" I*/
     //#include "petsc-private/kspimpl.h"
  #endif
#else
  #include "private/kspimpl.h"   /*I "petscksp.h" I*/
#endif

#include <StGermain/libStGermain/src/StGermain.h>
#include <StgDomain/libStgDomain/src/StgDomain.h>
#include <StgFEM/libStgFEM/src/StgFEM.h>
#include <PICellerator/libPICellerator/src/PICellerator.h>
#include <Underworld/libUnderworld/src/Underworld.h>
#include "Solvers/SLE/src/SLE.h" /* to give the AugLagStokes_SLE type */
#include "Solvers/KSPSolvers/src/KSPSolvers.h"

#include "BSSCR.h"
#include "createK2.h"
#include "writeMatVec.h"
#include "ksp_pressure_nullspace.h"

#define KSPBSSCR         "bsscr"

EXTERN_C_BEGIN
EXTERN PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_BSSCR(KSP);
EXTERN_C_END

PetscErrorCode BSSCR_DRIVER_auglag( KSP ksp, Mat stokes_A, Vec stokes_x, Vec stokes_b, Mat approxS,
                                    MatStokesBlockScaling BA, PetscTruth sym, KSP_BSSCR * bsscrp_self );
PetscErrorCode BSSCR_DRIVER_flex( KSP ksp, Mat stokes_A, Vec stokes_x, Vec stokes_b, Mat approxS, KSP ksp_K,
                                  MatStokesBlockScaling BA, PetscTruth sym, KSP_BSSCR * bsscrp_self );

/*************************************************************************************************/
/****** START Helper Functions *******************************************************************/
#define BSSCR_GetPetscMatrix( matrix ) ( (Mat)(matrix) )
/****** END Helper Functions *********************************************************************/
/*************************************************************************************************/
typedef struct {
    void *ctx;
    KSP_BSSCR * bsscr;
} BSSCR_KSPConverged_Ctx;

/*************************************************************************************************/
/*************************************************************************************************/

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPSetConvergenceMinIts"
PetscErrorCode BSSCR_KSPSetConvergenceMinIts(KSP ksp, PetscInt n, KSP_BSSCR * bsscr)
{
      BSSCR_KSPConverged_Ctx *ctx;
      PetscErrorCode ierr;

      PetscFunctionBegin;

      bsscr->min_it = n; /* set minimum its */

#if ( (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5 ) )
      ierr = Stg_PetscNew(BSSCR_KSPConverged_Ctx,&ctx);CHKERRQ(ierr);
      ierr = KSPConvergedDefaultCreate(&ctx->ctx);CHKERRQ(ierr);
      ctx->bsscr=bsscr;
      ierr = KSPSetConvergenceTest(ksp,BSSCR_KSPConverged,ctx,BSSCR_KSPConverged_Destroy);CHKERRQ(ierr);
#endif
#if ( (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR <=4 ) )
      ierr = Stg_PetscNew(BSSCR_KSPConverged_Ctx,&ctx);CHKERRQ(ierr);
      ierr = KSPDefaultConvergedCreate(&ctx->ctx);CHKERRQ(ierr);
      ctx->bsscr=bsscr;
      ierr = KSPSetConvergenceTest(ksp,BSSCR_KSPConverged,ctx,BSSCR_KSPConverged_Destroy);CHKERRQ(ierr);
#endif
#if ( PETSC_VERSION_MAJOR < 3)
      ierr = KSPSetConvergenceTest(ksp,BSSCR_KSPConverged,(void*)bsscr);CHKERRQ(ierr);
#endif

      PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPConverged"
PetscErrorCode BSSCR_KSPConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
{
  PetscErrorCode           ierr;
#if(PETSC_VERSION_MAJOR == 3)
  BSSCR_KSPConverged_Ctx  *ctx = (BSSCR_KSPConverged_Ctx *)cctx;
  KSP_BSSCR               *bsscr = ctx->bsscr;
#else
  KSP_BSSCR               *bsscr = (KSP_BSSCR*)cctx;
#endif

  PetscFunctionBegin;
#if ( (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >=5 ) )
  ierr = KSPConvergedDefault(ksp,n,rnorm,reason,ctx->ctx);CHKERRQ(ierr);
#endif
#if ( (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR <=4 ) )
  ierr = KSPDefaultConverged(ksp,n,rnorm,reason,ctx->ctx);CHKERRQ(ierr);
#endif
#if ( PETSC_VERSION_MAJOR < 3)
  ierr = KSPDefaultConverged(ksp,n,rnorm,reason,cctx);CHKERRQ(ierr);
#endif
  if (*reason) {
    ierr = PetscInfo2(ksp,"default convergence test KSP iterations=%D, rnorm=%G\n",n,rnorm);CHKERRQ(ierr);
  }
  if(ksp->its < bsscr->min_it){
      ksp->reason          = KSP_CONVERGED_ITERATING;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPConverged_Destroy"
PetscErrorCode BSSCR_KSPConverged_Destroy(void *cctx)
{
  BSSCR_KSPConverged_Ctx *ctx = (BSSCR_KSPConverged_Ctx *)cctx;
  PetscErrorCode           ierr;

  PetscFunctionBegin;
#if ( (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >=5 ) )
  ierr = KSPConvergedDefaultDestroy(ctx->ctx);CHKERRQ(ierr);
#else
  ierr = KSPDefaultConvergedDestroy(ctx->ctx);CHKERRQ(ierr);
#endif
  ierr = PetscFree(ctx);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPRegisterBSSCR"
PetscErrorCode PETSCKSP_DLLEXPORT KSPRegisterBSSCR(const char path[])
{
    PetscErrorCode ierr;

    PetscFunctionBegin;
    ierr = Stg_KSPRegister(KSPBSSCR, path, "KSPCreate_BSSCR", KSPCreate_BSSCR );CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPSolve_BSSCR"
PetscErrorCode  KSPSolve_BSSCR(KSP ksp)
{
    Mat            Amat,Pmat; /* Stokes Matrix and it's Preconditioner matrix: Both 2x2 PetscExt block matrices */
    Vec            B, X; /* rhs and solution vectors */
    //MatStructure   pflag;
    PetscErrorCode ierr;
    KSP_BSSCR *    bsscr;
    Stokes_SLE *    SLE;
    //PETScMGSolver * MG;
    Mat K,D,ApproxS;
    MatStokesBlockScaling BA;
    PetscTruth flg, sym, augment;
    double TotalSolveTime;

    PetscFunctionBegin;

    TotalSolveTime = MPI_Wtime();

    PetscPrintf( PETSC_COMM_WORLD, "\nBSSCR -- Block Stokes Schur Compliment Reduction Solver \n");
    /** Get the stokes Block matrix and its preconditioner matrix */
    ierr = Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,PETSC_NULL);CHKERRQ(ierr);

    /** In Petsc proper, KSP's ksp->data is usually set in KSPCreate_XXX function.
        Here it is set in the _StokesBlockKSPInterface_Solve function instead so that we can ensure that the solver
        has everything it needs */

    bsscr         = (KSP_BSSCR*)ksp->data;
    //MG            = (PETScMGSolver*)bsscr->mg;
    SLE           = (Stokes_SLE*)bsscr->st_sle;
    X             = ksp->vec_sol;
    B             = ksp->vec_rhs;

    if( bsscr->do_scaling ){
        (*bsscr->scale)(ksp); /* scales everything including the UW preconditioner */
        BA =  bsscr->BA;
    }

    if( (bsscr->k2type != 0) ){
      if(bsscr->buildK2 != PETSC_NULL) {
            (*bsscr->buildK2)(ksp); /* building K2 from scaled version of stokes operators: K2 lives on bsscr struct = ksp->data */
        }
    }

    /* get sub matrix / vector objects */
    MatNestGetSubMat( Amat, 0,0, &K );
    /* Underworld preconditioner matrix*/
    ApproxS = PETSC_NULL;
    if( ((StokesBlockKSPInterface*)SLE->solver)->preconditioner ) { /* SLE->solver->st_sle == SLE here, by the way */
        StiffnessMatrix *preconditioner;
        preconditioner = ((StokesBlockKSPInterface*)SLE->solver)->preconditioner;
        ApproxS = BSSCR_GetPetscMatrix( preconditioner->matrix );
    }

    sym = bsscr->DIsSym;

    MatNestGetSubMat( Amat, 1,0, &D );if(!D){ PetscPrintf( PETSC_COMM_WORLD, "D does not exist but should!!\n"); exit(1); }

    /**********************************************************/
    /******* SOLVE!! ******************************************/
    /**********************************************************/

    flg = PETSC_FALSE;
    augment = PETSC_TRUE;
    PetscOptionsGetTruth(PETSC_NULL, "-augmented_lagrangian", &augment, &flg);
    BSSCR_DRIVER_auglag( ksp, Amat, X, B, ApproxS, BA, sym, bsscr );

    /**********************************************************/
    /***** END SOLVE!! ****************************************/
    /**********************************************************/
    if( bsscr->do_scaling ){
        (*bsscr->unscale)(ksp);  }
    if( (bsscr->k2type != 0) && bsscr->K2 != PETSC_NULL ){
        if(bsscr->k2type != K2_SLE){/* don't destroy here, as in this case, K2 is just pointing to an existing matrix on the SLE */
            Stg_MatDestroy(&bsscr->K2 );
        }
        bsscr->K2built = PETSC_FALSE;
        bsscr->K2 = PETSC_NULL;  }

    ksp->reason = KSP_CONVERGED_RTOL;

    TotalSolveTime =  MPI_Wtime() - TotalSolveTime;
    PetscPrintf( PETSC_COMM_WORLD, "  Total BSSCR Linear solve time: %lf seconds\n\n", TotalSolveTime);
    bsscr->solver->stats.total_time=TotalSolveTime;
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPDestroy_BSSCR"
PetscErrorCode KSPDestroy_BSSCR(KSP ksp)
{
    KSP_BSSCR     *bsscr = (KSP_BSSCR *)ksp->data;
    Mat            K2=bsscr->K2;
    Vec            t,v;
    MatStokesBlockScaling BA=bsscr->BA;
    PetscErrorCode ierr;

    PetscFunctionBegin;

    t=bsscr->t;
    v=bsscr->v;
    if ( v ){ Stg_VecDestroy(&v); }
    if ( t ){ Stg_VecDestroy(&t); }
    if( K2 ){ Stg_MatDestroy(&K2 ); }/* shouldn't need this now */
    if(BA) BSSCR_MatStokesBlockScalingDestroy( BA );
    ierr = PetscFree(ksp->data);CHKERRQ(ierr);

    PetscFunctionReturn(0);
}
static const char *K2Types[] = {"NULL","DGMGD","GMG","GG","SLE","K2Type","K2_",0};
#undef __FUNCT__
#define __FUNCT__ "KSPSetFromOptions_BSSCR"
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR <6) )
PetscErrorCode KSPSetFromOptions_BSSCR(KSP ksp)
{

    PetscTruth  flg;
    KSP_BSSCR  *bsscr = (KSP_BSSCR *)ksp->data;
    PetscErrorCode ierr;
    PetscFunctionBegin;
    ierr = PetscOptionsHead("KSP BSSCR options");CHKERRQ(ierr);
    /* if this ksp has a prefix "XXX_" it will be automatically added to the options. e.g. -ksp_test -> -XXX_ksp_test */
    /* ierr = PetscOptionsTruth("-ksp_test","Test KSP flag","nil",PETSC_FALSE,&test,PETSC_NULL);CHKERRQ(ierr); */
    /* if(test){ PetscPrintf( PETSC_COMM_WORLD,  "\n\n-----  test flag set  ------\n\n"); } */
    ierr = PetscOptionsEnum("-ksp_k2_type","Augmented Lagrangian matrix type","",K2Types, bsscr->k2type,(PetscEnum*)&bsscr->k2type,&flg);CHKERRQ(ierr);
    //if(flg){  PetscPrintf( PETSC_COMM_WORLD,  "-----  k2 type is  ------\n"); }
    ierr = PetscOptionsTail();CHKERRQ(ierr);
    PetscFunctionReturn(0);
}
#else
PetscErrorCode KSPSetFromOptions_BSSCR(Stg_PetscOptions *PetscOptionsObject, KSP ksp)
{

    PetscTruth  flg;
    KSP_BSSCR  *bsscr = (KSP_BSSCR *)ksp->data;
    PetscErrorCode ierr;
    PetscFunctionBegin;
    ierr = PetscOptionsHead(PetscOptionsObject,"KSP BSSCR options");CHKERRQ(ierr);
    /* if this ksp has a prefix "XXX_" it will be automatically added to the options. e.g. -ksp_test -> -XXX_ksp_test */
    /* ierr = PetscOptionsTruth("-ksp_test","Test KSP flag","nil",PETSC_FALSE,&test,PETSC_NULL);CHKERRQ(ierr); */
    /* if(test){ PetscPrintf( PETSC_COMM_WORLD,  "\n\n-----  test flag set  ------\n\n"); } */
    ierr = PetscOptionsEnum("-ksp_k2_type","Augmented Lagrangian matrix type","",K2Types, bsscr->k2type,(PetscEnum*)&bsscr->k2type,&flg);CHKERRQ(ierr);
    //if(flg){  PetscPrintf( PETSC_COMM_WORLD,  "-----  k2 type is  ------\n"); }
    ierr = PetscOptionsTail();CHKERRQ(ierr);
    PetscFunctionReturn(0);
}
#endif

#undef __FUNCT__
#define __FUNCT__ "KSPView_BSSCR"
PetscErrorCode KSPView_BSSCR(KSP ksp,PetscViewer viewer)
{
    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPSetUp_BSSCR"
PetscErrorCode KSPSetUp_BSSCR(KSP ksp)
{
    KSP_BSSCR  *bsscr = (KSP_BSSCR *)ksp->data;
    Mat K;
    Stokes_SLE*  stokesSLE  = (Stokes_SLE*)bsscr->st_sle;
    PetscTruth ismumps,augment,scale,konly,found,conp,checkerp;

    PetscFunctionBegin;

    BSSCR_PetscExtStokesSolversInitialize();

    K = stokesSLE->kStiffMat->matrix;
    bsscr->K2=PETSC_NULL;

    BSSCR_MatStokesBlockScalingCreate( &(bsscr->BA) );/* allocate memory for scaling struct */
    found = PETSC_FALSE;
    augment = PETSC_TRUE;
    PetscOptionsGetTruth(PETSC_NULL, "-augmented_lagrangian", &augment, &found);
    if(augment){
        bsscr->buildK2 = bsscr_buildK2; }
    else {
        bsscr->buildK2 = PETSC_NULL; }

    /***************************************************************************************************************/
    /** Do scaling *************************************************************************************************/
    /***************************************************************************************************************/
    found = PETSC_FALSE;
    scale = PETSC_FALSE;/* scaling off by default */
    PetscOptionsGetTruth(PETSC_NULL, "-rescale_equations", &scale, &found);
    Stg_PetscObjectTypeCompare((PetscObject)K, "mataijmumps", &ismumps);/** older versions of petsc have this */
    if(ismumps && scale){
        PetscPrintf( PETSC_COMM_WORLD, "\t* Not applying scaling to matrices as MatGetRowMax operation not defined for MATAIJMUMPS matrix type \n");
        scale = PETSC_FALSE;
    }
    if( scale ) {
        bsscr->scale       = KSPScale_BSSCR;
        bsscr->unscale     = KSPUnscale_BSSCR;
        bsscr->do_scaling  = PETSC_TRUE;
        bsscr->scaletype   = KONLY;
        konly = PETSC_TRUE;
        found = PETSC_FALSE;
        PetscOptionsGetTruth(PETSC_NULL, "-k_scale_only", &konly, &found);
        if(!konly){ bsscr->scaletype = DEFAULT; }
    }

    /***************************************************************************************************************/
    /**  Set up functions for building Pressure Null Space Vectors *************************************************/
    /***************************************************************************************************************/
    found = PETSC_FALSE;
    checkerp = PETSC_FALSE;
    PetscOptionsGetTruth(PETSC_NULL, "-remove_checkerboard_pressure_null_space", &checkerp, &found );
    if(checkerp){
        bsscr->check_cb_pressureNS    = PETSC_TRUE;
        bsscr->check_pressureNS       = PETSC_TRUE;
        bsscr->buildPNS               = KSPBuildPressure_CB_Nullspace_BSSCR;
    }
    found = PETSC_FALSE;
    conp = PETSC_FALSE;
    PetscOptionsGetTruth(PETSC_NULL, "-remove_constant_pressure_null_space", &conp, &found );
    if(conp){
        bsscr->check_const_pressureNS = PETSC_TRUE;
        bsscr->check_pressureNS       = PETSC_TRUE;
        bsscr->buildPNS               = KSPBuildPressure_Const_Nullspace_BSSCR;
    }
    /***************************************************************************************************************/
    PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "KSPCreate_BSSCR"
PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_BSSCR(KSP ksp)
{
    KSP_BSSCR  *bsscr;
    PetscErrorCode ierr;
    PetscFunctionBegin;
    ierr = Stg_PetscNew(KSP_BSSCR,&bsscr);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)ksp,sizeof(KSP_BSSCR));CHKERRQ(ierr);
    //ierr = PetscNewLog(ksp,KSP_BSSCR,&bsscr);CHKERRQ(ierr);
    ksp->data                              = (void*)bsscr;

    #if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >= 2 ) )
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,0);CHKERRQ(ierr);
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);CHKERRQ(ierr);
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,0);CHKERRQ(ierr);
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
    #endif
    /*
       Sets the functions that are associated with this data structure
       (in C++ this is the same as defining virtual functions)
    */
    ksp->ops->setup                = KSPSetUp_BSSCR;
    ksp->ops->solve                = KSPSolve_BSSCR;
    ksp->ops->destroy              = KSPDestroy_BSSCR;
    ksp->ops->view                 = KSPView_BSSCR;
    ksp->ops->setfromoptions       = KSPSetFromOptions_BSSCR;
    ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    ksp->ops->buildresidual        = KSPDefaultBuildResidual;

//    bsscr->k2type=K2_GMG;
    bsscr->k2type      = 0;
    bsscr->do_scaling  = PETSC_FALSE;
    bsscr->scaled      = PETSC_FALSE;
    bsscr->K2built     = PETSC_FALSE;
    bsscr->check_cb_pressureNS     = PETSC_FALSE;/* checker board nullspace */
    bsscr->check_const_pressureNS  = PETSC_FALSE;/* constant nullspace */
    bsscr->check_pressureNS        = PETSC_FALSE;/* is true if either of above two are true */
    bsscr->t           = NULL;/* null space vectors for pressure */
    bsscr->v           = NULL;
    bsscr->nstol       = 1e-7;/* null space detection tolerance */
    bsscr->uStar       = NULL;
    bsscr->been_here   = 0;
    PetscFunctionReturn(0);
}
EXTERN_C_END


text:
static char help[] =
"Solves 2D advection-diffusion problems using FD discretization,\n"
"structured-grid (DMDA), and -snes_fd_color.  Option prefix -bth_.\n"
"Equation:\n"
"    - eps Laplacian u + Div (a(x,y) u) = g(x,y),\n"
"where the (vector) wind a(x,y) and (scalar) source g(x,y) are given smooth\n"
"functions.  The domain is S = (-1,1)^2 with Dirichlet boundary conditions:\n"
"    u = b(x,y) on boundary S\n"
"where b(x,y) is a given smooth function.  Problems include: NOWIND, LAYER,\n"
"and GLAZE.  The first of these has a=0 while LAYER and GLAZE are\n"
"Examples 6.1.1 and 6.1.4 in Elman et al (2014), respectively.\n"
"Advection can be discretized by first-order upwinding (none), centered, or a\n"
"van Leer limiter scheme.  Option allows switching to none limiter on all grids\n"
"for which the mesh Peclet P^h exceeds a threshold (default: 1).\n\n";

#include <petsc.h>

typedef enum {NONE, CENTERED, VANLEER} LimiterType;
static const char *LimiterTypes[] = {"none","centered","vanleer",
                                     "LimiterType", "", NULL};

static PetscReal centered(PetscReal theta) {
    return 0.5;
}

static PetscReal vanleer(PetscReal theta) {
    const PetscReal abstheta = PetscAbsReal(theta);
    return 0.5 * (theta + abstheta) / (1.0 + abstheta);   // 4 flops
}

typedef PetscReal (*LimiterFcn)(PetscReal);

static LimiterFcn limiterptr[] = {NULL, &centered, &vanleer};

typedef enum {NOWIND, LAYER, GLAZE} ProblemType;
static const char *ProblemTypes[] = {"nowind", "layer", "glaze",
                                     "ProblemType", "", NULL};

typedef struct {
    ProblemType  problem;
    PetscReal    eps,                              // diffusion eps > 0
                 a_scale,                          // scale for wind
                 peclet_threshold,
                 (*limiter_fcn)(PetscReal),
                 (*g_fcn)(PetscReal, PetscReal, void*),  // right-hand-side source
                 (*b_fcn)(PetscReal, PetscReal, void*);  // boundary condition
    PetscBool    none_on_peclet,                   // if true use none limiter when P^h > threshold
                 small_peclet_achieved;            // true if on finest grid P^h <= threshold
} AdCtx;

// used for source functions
static PetscReal zero(PetscReal x, PetscReal y, void *user) {
    return 0.0;
}

// problem NOWIND: same problem as ./fish -fsh_problem manuexp
static PetscReal nowind_u(PetscReal x, PetscReal y, void *user) {  // exact solution
    return - x * PetscExpReal(y);
}

static PetscReal nowind_g(PetscReal x, PetscReal y, void *user) {
    AdCtx* usr = (AdCtx*)user;
    return usr->eps * x * PetscExpReal(y);
}

static PetscReal nowind_b(PetscReal x, PetscReal y, void *user) {
    return nowind_u(x,y,user);
}

// problem LAYER:  Elman page 237, Example 6.1.1
static PetscReal layer_u(PetscReal x, PetscReal y, void *user) {  // exact solution
    AdCtx* usr = (AdCtx*)user;
    return x * (1.0 - PetscExpReal((y-1) / usr->eps))
             / (1.0 - PetscExpReal(- 2.0 / usr->eps));
}

static PetscReal layer_b(PetscReal x, PetscReal y, void *user) {
    return layer_u(x,y,user);
}

// problem GLAZE:  Elman page 240, Example 6.1.4
static PetscReal glaze_b(PetscReal x, PetscReal y, void *user) {
    if (x > 0.0 && y < x && y > -x)
       return 1.0;   // along x=1 boundary
    else
       return 0.0;
}

typedef PetscReal (*PointwiseFcn)(PetscReal,PetscReal,void*);

static PointwiseFcn uexptr[] = {&nowind_u, &layer_u, NULL};
static PointwiseFcn gptr[]   = {&nowind_g, &zero,    &zero};
static PointwiseFcn bptr[]   = {&nowind_b, &layer_b, &glaze_b};

/* This vector function returns q=0,1 component.  It is used in
   FormFunctionLocal() to get a(x,y). */
static PetscReal wind_a(PetscReal x, PetscReal y, PetscInt q, AdCtx *user) {
    switch (user->problem) {
        case NOWIND:
            return 0.0;
        case LAYER:
            return (q == 0) ? 0.0 : 1.0;
        case GLAZE:
            return (q == 0) ? 2.0*y*(1.0-x*x) : -2.0*x*(1.0-y*y);
        default:
            return NAN;
    }
}

extern PetscErrorCode FormUExact(DMDALocalInfo*,AdCtx*,
                                 PetscReal (*)(PetscReal, PetscReal, void*),Vec);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscReal**,PetscReal**,AdCtx*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da, da_after;
    SNES           snes;
    Vec            u_initial, u;
    DMDALocalInfo  info;
    PointwiseFcn   uexact_fcn;
    LimiterType    limiter = NONE;
    PetscBool      init_exact = PETSC_FALSE;
    AdCtx          user;

    ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

    user.eps = 0.005;
    user.none_on_peclet = PETSC_FALSE;
    user.small_peclet_achieved = PETSC_FALSE;
    user.problem = LAYER;
    user.a_scale = 1.0;   // this could be made dependent on problem
    user.peclet_threshold = 1.0;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"bth_",
               "both (2D advection-diffusion solver) options",""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-eps","positive diffusion coefficient",
               "both.c",user.eps,&(user.eps),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-init_exact","use exact solution for initialization",
               "both.c",init_exact,&init_exact,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-limiter","flux-limiter type",
               "both.c",LimiterTypes,
               (PetscEnum)limiter,(PetscEnum*)&limiter,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-none_on_peclet",
               "on coarse grids such that mesh Peclet P^h exceeds threshold, switch to none limiter",
               "both.c",user.none_on_peclet,&(user.none_on_peclet),NULL);
               CHKERRQ(ierr);
    ierr = PetscOptionsReal("-peclet_threshold",
               "if mesh Peclet P^h is above this value, switch to none (used with -bth_none_on_peclet)",
               "both.c",user.peclet_threshold,&(user.peclet_threshold),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-problem","problem type",
               "both.c",ProblemTypes,
               (PetscEnum)(user.problem),(PetscEnum*)&(user.problem),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    if (user.eps <= 0.0) {
        SETERRQ1(PETSC_COMM_SELF,1,"eps=%.3f invalid ... eps > 0 required",user.eps);
    }
    user.limiter_fcn = limiterptr[limiter];
    uexact_fcn = uexptr[user.problem];
    user.g_fcn = gptr[user.problem];
    user.b_fcn = bptr[user.problem];

    ierr = DMDACreate2d(PETSC_COMM_WORLD,
        DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
        3,3,                          // default to hx=hy=1 grid
        PETSC_DECIDE,PETSC_DECIDE,
        1,2,                          // d.o.f, stencil width
        NULL,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    if (user.problem == NOWIND) {
        ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,-1.0,1.0); CHKERRQ(ierr);
    } else {
        ierr = DMDASetUniformCoordinates(da,-1.0,1.0,-1.0,1.0,-1.0,1.0); CHKERRQ(ierr);
    }
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
    ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
            (DMDASNESFunction)FormFunctionLocal,&user);CHKERRQ(ierr);
    ierr = SNESSetApplicationContext(snes,&user); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

    ierr = DMGetGlobalVector(da,&u_initial); CHKERRQ(ierr);
    if (init_exact) {
        ierr = FormUExact(&info,&user,uexact_fcn,u_initial); CHKERRQ(ierr);
    } else {
        ierr = VecSet(u_initial,0.0); CHKERRQ(ierr);
    }
    ierr = SNESSolve(snes,NULL,u_initial); CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(da,&u_initial); CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);

    ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
    if (user.none_on_peclet && !user.small_peclet_achieved) {
         ierr = PetscPrintf(PETSC_COMM_WORLD,
             "WARNING: -bth_none_on_peclet set but finest grid used NONE limiter\n");
             CHKERRQ(ierr);
    }
    ierr = PetscPrintf(PETSC_COMM_WORLD,
         "done on %d x %d grid (problem = %s, eps = %g, limiter = %s)\n",
         info.mx,info.my,ProblemTypes[user.problem],user.eps,LimiterTypes[limiter]);
         CHKERRQ(ierr);

    if (uexact_fcn != NULL) {
        Vec     u_exact;
        PetscReal  xymin[2], xymax[2], hx, hy, err2;
        ierr = DMCreateGlobalVector(da_after,&u_exact); CHKERRQ(ierr);
        ierr = FormUExact(&info,&user,uexact_fcn,u_exact); CHKERRQ(ierr);
        ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr);    // u <- u + (-1.0) u_exact
        ierr = VecNorm(u,NORM_2,&err2); CHKERRQ(ierr);
        ierr = DMGetBoundingBox(da_after,xymin,xymax); CHKERRQ(ierr);
        hx = (xymax[0] - xymin[0]) / (info.mx - 1);
        hy = (xymax[1] - xymin[1]) / (info.my - 1);
        err2 *= PetscSqrtReal(hx * hy);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
             "numerical error:  |u-uexact|_2 = %.4e\n",err2); CHKERRQ(ierr);
        ierr = VecDestroy(&u_exact); CHKERRQ(ierr);
    }

    ierr = SNESDestroy(&snes); CHKERRQ(ierr);
    return PetscFinalize();
}

PetscErrorCode FormUExact(DMDALocalInfo *info, AdCtx *usr,
                          PetscReal (*uexact)(PetscReal, PetscReal, void*), Vec uex) {
    PetscErrorCode  ierr;
    PetscInt        i, j;
    PetscReal       xymin[2], xymax[2], hx, hy, x, y, **auex;

    if (uexact == NULL) {
        SETERRQ(PETSC_COMM_SELF,1,"exact solution not available");
    }
    ierr = DMGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    ierr = DMDAVecGetArray(info->da, uex, &auex);CHKERRQ(ierr);
    for (j=info->ys; j<info->ys+info->ym; j++) {
        y = xymin[1] + j * hy;
        for (i=info->xs; i<info->xs+info->xm; i++) {
            x = xymin[0] + i * hx;
            auex[j][i] = (*uexact)(x,y,usr);
        }
    }
    ierr = DMDAVecRestoreArray(info->da, uex, &auex);CHKERRQ(ierr);
    return 0;
}

/* compute residuals:
     F_ij = hx * hy * (- eps Laplacian u + Div (a(x,y) u) - g(x,y))
at boundary points:
     F_ij = c (u - b(x,y))
note compass notation for faces of cells; advection scheme evaluates flux at
each E,N face once and then includes it in two residuals:
     N
  -------
  |     |
W |  *  | E
  |     |
  -------
     S
*/
PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscReal **au,
                                 PetscReal **aF, AdCtx *usr) {
    PetscErrorCode ierr;
    PetscInt        i, j, p;
    PetscReal       xymin[2], xymax[2], hx, hy, Ph, hx2, hy2, scF, scBC,
                    x, y, uE, uW, uN, uS, uxx, uyy,
                    ap, flux, u_up, u_dn, u_far, theta;
    PetscReal       (*limiter)(PetscReal);
    PetscBool       iowned, jowned, ip1owned, jp1owned;
    PetscLogDouble  ff;

    ierr = DMGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    limiter = usr->limiter_fcn;
    if (usr->none_on_peclet) {
        Ph = usr->a_scale * PetscMax(hx,hy) / usr->eps;  // mesh Peclet number
        if (Ph > usr->peclet_threshold)
            limiter = NULL;
        else
            usr->small_peclet_achieved = PETSC_TRUE;
    }
    hx2 = hx * hx;
    hy2 = hy * hy;
    scF = hx * hy;  // scale residuals
    scBC = scF * usr->eps * 2.0 * (1.0 / hx2 + 1.0 / hy2); // scale b.c. residuals

    // for owned cells, compute non-advective parts of residual at cell center
    for (j=info->ys; j<info->ys+info->ym; j++) {
        y = xymin[1] + j * hy;
        for (i=info->xs; i<info->xs+info->xm; i++) {
            x = xymin[0] + i * hx;
            if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {
                aF[j][i] = scBC * (au[j][i] - (*usr->b_fcn)(x,y,usr));
            } else {
                uE = (i+1 == info->mx-1) ? (*usr->b_fcn)(xymax[0],y,usr) : au[j][i+1];
                uW = (i-1 == 0)          ? (*usr->b_fcn)(xymin[0],y,usr) : au[j][i-1];
                uxx = (uE - 2.0 * au[j][i] + uW) / hx2;
                uN = (j+1 == info->my-1) ? (*usr->b_fcn)(x,xymax[1],usr) : au[j+1][i];
                uS = (j-1 == 0)          ? (*usr->b_fcn)(x,xymin[1],usr) : au[j-1][i];
                uyy = (uN - 2.0 * au[j][i] + uS) / hy2;
                aF[j][i] = scF * (- usr->eps * (uxx + uyy) - (*usr->g_fcn)(x,y,usr));
            }
        }
    }
    ierr = PetscLogFlops(14.0*info->xm*info->ym); CHKERRQ(ierr);

    // for each E,N face of an *owned* cell at (x,y) and (i,j), compute flux at
    //     the face center and then add that to the correct residual
    // note start offset of -1; gets W,S faces of owned cells living on ownership
    //     boundaries for i,j resp.
    // there are (xm+1)*(ym+1)*2 fluxes to evaluate
    for (j=info->ys-1; j<info->ys+info->ym; j++) {
        y = xymin[1] + j * hy;
        // if y<0 or y=1 at cell center then no need to compute *any* E,N face-center fluxes
        if (j < 0 || j == info->my-1)
            continue;
        for (i=info->xs-1; i<info->xs+info->xm; i++) {
            x = xymin[0] + i * hx;
            // if x<0 or x=1 at cell center then no need to compute *any* E,N face-center fluxes
            if (i < 0 || i == info->mx-1)
                continue;
            // get E (p=0) and N (p=1) cell face-center flux contributions
            for (p = 0; p < 2; p++) {
                // get pth component of wind and locations determined by wind direction
                ap = (p == 0) ? wind_a(x+hx/2.0,y,p,usr) : wind_a(x,y+hy/2.0,p,usr);
                if (p == 0)
                    if (ap >= 0.0) {
                        u_up  = (i == 0)            ? (*usr->b_fcn)(xymin[0],y,usr) : au[j][i];
                        u_dn  = (i+1 == info->mx-1) ? (*usr->b_fcn)(xymax[0],y,usr) : au[j][i+1];
                        u_far = (i-1 <= 0)          ? (*usr->b_fcn)(xymin[0],y,usr) : au[j][i-1];
                    } else {
                        u_up  = (i+1 == info->mx-1) ? (*usr->b_fcn)(xymax[0],y,usr) : au[j][i+1];
                        u_dn  = (i == 0)            ? (*usr->b_fcn)(xymin[0],y,usr) : au[j][i];
                        u_far = (i+2 >= info->mx-1) ? (*usr->b_fcn)(xymax[0],y,usr) : au[j][i+2];
                    }
                else  // p == 1
                    if (ap >= 0.0) {
                        u_up  = (j == 0)            ? (*usr->b_fcn)(x,xymin[1],usr) : au[j][i];
                        u_dn  = (j+1 == info->my-1) ? (*usr->b_fcn)(x,xymax[1],usr) : au[j+1][i];
                        u_far = (j-1 <= 0)          ? (*usr->b_fcn)(x,xymin[1],usr) : au[j-1][i];
                    } else {
                        u_up  = (j+1 == info->my-1) ? (*usr->b_fcn)(x,xymax[1],usr) : au[j+1][i];
                        u_dn  = (j == 0)            ? (*usr->b_fcn)(x,xymin[1],usr) : au[j][i];
                        u_far = (j+2 >= info->my-1) ? (*usr->b_fcn)(x,xymax[1],usr) : au[j+2][i];
                    }
                // first-order upwind flux plus correction if have limiter
                flux = ap * u_up;
                if (limiter != NULL && u_dn != u_up) {
                    theta = (u_up - u_far) / (u_dn - u_up);
                    flux += ap * (*limiter)(theta) * (u_dn - u_up);
                }
                // update non-boundary and owned residual F_ij on both sides of computed flux
                // note: 1) aF[] does not have stencil width, 2) F_ij is scaled by scF = hx * hy
                iowned = (i >= info->xs);
                jowned = (j >= info->ys);
                ip1owned = (i+1 < info->xs + info->xm);
                jp1owned = (j+1 < info->ys + info->ym);
                switch (p) {
                    case 0:  // flux at E
                        if (i > 0 && j > 0 && iowned && jowned)
                            aF[j][i] += hy * flux;  // flux out of i,j at E
                        if (i+1 < info->mx-1 && j > 0 && ip1owned && jowned)
                            aF[j][i+1] -= hy * flux;  // flux into i+1,j at W
                        break;
                    case 1:  // flux at N
                        if (j > 0 && i > 0 && iowned && jowned)
                            aF[j][i] += hx * flux;  // flux out of i,j at N
                        if (j+1 < info->my-1 && i > 0 && iowned && jp1owned)
                            aF[j+1][i] -= hx * flux;  // flux into i,j+1 at S
                        break;
                }
            }
        }
    }
    // ff = flops per flux evaluation
    ff = (limiter == NULL) ? 6.0 : 13.0;
    if (limiter == &vanleer)
        ff += 4.0;
    ierr = PetscLogFlops(ff*2.0*(1.0+info->xm)*(1.0+info->ym)); CHKERRQ(ierr);
    return 0;
}