text:
/* linalg/gsl_linalg.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007, 2019 Gerard Jungman, Brian Gough, Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_LINALG_H__
#define __GSL_LINALG_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_blas.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
#define __BEGIN_DECLS extern "C" {
#define __END_DECLS }
#else
#define __BEGIN_DECLS           /* empty */
#define __END_DECLS             /* empty */
#endif

__BEGIN_DECLS

typedef enum
  {
    GSL_LINALG_MOD_NONE = 0,
    GSL_LINALG_MOD_TRANSPOSE = 1,
    GSL_LINALG_MOD_CONJUGATE = 2
  }
gsl_linalg_matrix_mod_t;

/* Note: You can now use the gsl_blas_dgemm function instead of matmult */

/* Simple implementation of matrix multiply.
 * Calculates C = A.B
 *
 * exceptions: GSL_EBADLEN
 */
GSL_FUN int gsl_linalg_matmult (const gsl_matrix * A,
                        const gsl_matrix * B,
                        gsl_matrix * C);

/* Simple implementation of matrix multiply.
 * Allows transposition of either matrix, so it
 * can compute A.B or Trans(A).B or A.Trans(B) or Trans(A).Trans(B)
 *
 * exceptions: GSL_EBADLEN
 */
GSL_FUN int gsl_linalg_matmult_mod (const gsl_matrix * A,
                            gsl_linalg_matrix_mod_t modA,
                            const gsl_matrix * B,
                            gsl_linalg_matrix_mod_t modB,
                            gsl_matrix * C);

/* Calculate the matrix exponential by the scaling and
 * squaring method described in Moler + Van Loan,
 * SIAM Rev 20, 801 (1978). The mode argument allows
 * choosing an optimal strategy, from the table
 * given in the paper, for a given precision.
 *
 * exceptions: GSL_ENOTSQR, GSL_EBADLEN
 */
GSL_FUN int gsl_linalg_exponential_ss(
  const gsl_matrix * A,
  gsl_matrix * eA,
  gsl_mode_t mode
  );

/* Householder Transformations */

GSL_FUN double gsl_linalg_householder_transform (gsl_vector * v);
GSL_FUN double gsl_linalg_householder_transform2 (double * alpha, gsl_vector * v);
GSL_FUN gsl_complex gsl_linalg_complex_householder_transform (gsl_vector_complex * v);

GSL_FUN int gsl_linalg_householder_hm (double tau,
                               const gsl_vector * v,
                               gsl_matrix * A);

GSL_FUN int gsl_linalg_householder_mh (double tau,
                               const gsl_vector * v,
                               gsl_matrix * A);

GSL_FUN int gsl_linalg_householder_hv (double tau,
                               const gsl_vector * v,
                               gsl_vector * w);

GSL_FUN int gsl_linalg_householder_left(const double tau,
                                const gsl_vector * v,
                                gsl_matrix * A,
                                gsl_vector * work);

GSL_FUN int gsl_linalg_householder_right(const double tau,
                                 const gsl_vector * v,
                                 gsl_matrix * A,
                                 gsl_vector * work);

GSL_FUN int gsl_linalg_householder_hm1 (double tau,
                                gsl_matrix * A);

GSL_FUN int gsl_linalg_complex_householder_hm (gsl_complex tau,
                                       const gsl_vector_complex * v,
                                       gsl_matrix_complex * A);

GSL_FUN int gsl_linalg_complex_householder_mh (gsl_complex tau,
                                       const gsl_vector_complex * v,
                                       gsl_matrix_complex * A);

GSL_FUN int gsl_linalg_complex_householder_hv (gsl_complex tau,
                                       const gsl_vector_complex * v,
                                       gsl_vector_complex * w);

GSL_FUN int gsl_linalg_complex_householder_left (const gsl_complex tau,
                                         const gsl_vector_complex * v,
                                         gsl_matrix_complex * A,
                                         gsl_vector_complex * work);

/* Hessenberg reduction */

GSL_FUN int gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau);
GSL_FUN int gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
                                 gsl_matrix * U);
GSL_FUN int gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
                                       gsl_matrix * U);
GSL_FUN int gsl_linalg_hessenberg_set_zero(gsl_matrix * H);
GSL_FUN int gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A,
                                    size_t top, gsl_vector *tau);

/* Hessenberg-Triangular reduction */

GSL_FUN int gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B,
                              gsl_matrix * U, gsl_matrix * V,
                              gsl_vector * work);

/* Singular Value Decomposition

 * exceptions:
 */

GSL_FUN int
gsl_linalg_SV_decomp (gsl_matrix * A,
                      gsl_matrix * V,
                      gsl_vector * S,
                      gsl_vector * work);

GSL_FUN int
gsl_linalg_SV_decomp_mod (gsl_matrix * A,
                          gsl_matrix * X,
                          gsl_matrix * V,
                          gsl_vector * S,
                          gsl_vector * work);

GSL_FUN int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A,
                                 gsl_matrix * Q,
                                 gsl_vector * S);

GSL_FUN int
gsl_linalg_SV_solve (const gsl_matrix * U,
                     const gsl_matrix * Q,
                     const gsl_vector * S,
                     const gsl_vector * b,
                     gsl_vector * x);

GSL_FUN int gsl_linalg_SV_leverage(const gsl_matrix *U, gsl_vector *h);

/* LU Decomposition, Gaussian elimination with partial pivoting
 */

GSL_FUN int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum);

GSL_FUN int gsl_linalg_LU_solve (const gsl_matrix * LU,
                         const gsl_permutation * p,
                         const gsl_vector * b,
                         gsl_vector * x);

GSL_FUN int gsl_linalg_LU_svx (const gsl_matrix * LU,
                       const gsl_permutation * p,
                       gsl_vector * x);

GSL_FUN int gsl_linalg_LU_refine (const gsl_matrix * A,
                          const gsl_matrix * LU,
                          const gsl_permutation * p,
                          const gsl_vector * b,
                          gsl_vector * x,
                          gsl_vector * work);

GSL_FUN int gsl_linalg_LU_invert (const gsl_matrix * LU,
                          const gsl_permutation * p,
                          gsl_matrix * inverse);
GSL_FUN int gsl_linalg_LU_invx (gsl_matrix * LU, const gsl_permutation * p);

GSL_FUN double gsl_linalg_LU_det (gsl_matrix * LU, int signum);
GSL_FUN double gsl_linalg_LU_lndet (gsl_matrix * LU);
GSL_FUN int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum);

/* Complex LU Decomposition */

GSL_FUN int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A,
                                  gsl_permutation * p,
                                  int *signum);

GSL_FUN int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU,
                                 const gsl_permutation * p,
                                 const gsl_vector_complex * b,
                                 gsl_vector_complex * x);

GSL_FUN int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU,
                               const gsl_permutation * p,
                               gsl_vector_complex * x);

GSL_FUN int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A,
                                  const gsl_matrix_complex * LU,
                                  const gsl_permutation * p,
                                  const gsl_vector_complex * b,
                                  gsl_vector_complex * x,
                                  gsl_vector_complex * work);

GSL_FUN int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU,
                                  const gsl_permutation * p,
                                  gsl_matrix_complex * inverse);
GSL_FUN int gsl_linalg_complex_LU_invx (gsl_matrix_complex * LU, const gsl_permutation * p);

GSL_FUN gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU,
                                       int signum);

GSL_FUN double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU);

GSL_FUN gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU,
                                          int signum);

/* QR decomposition */

GSL_FUN int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau);

GSL_FUN int gsl_linalg_QR_decomp_r (gsl_matrix * A, gsl_matrix * T);

GSL_FUN int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_solve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b,
                           gsl_vector * x, gsl_vector * residual);

GSL_FUN int gsl_linalg_QR_lssolve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b,
                             gsl_vector * x, gsl_vector * work);

GSL_FUN int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R, gsl_vector * w, const gsl_vector * v);

GSL_FUN int gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);

GSL_FUN int gsl_linalg_QR_QTvec_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);

GSL_FUN int gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);

GSL_FUN int gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);

GSL_FUN int gsl_linalg_QR_QTmat_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * B, gsl_matrix * work);

GSL_FUN int gsl_linalg_QR_matQ (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);

GSL_FUN int gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R);

GSL_FUN int gsl_linalg_QR_unpack_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q, gsl_matrix * R);

GSL_FUN int gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x);

GSL_FUN int gsl_linalg_QR_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work);

/* Q R P^T decomposition */

GSL_FUN int gsl_linalg_QRPT_decomp (gsl_matrix * A,
                            gsl_vector * tau,
                            gsl_permutation * p,
                            int *signum,
                            gsl_vector * norm);

GSL_FUN int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
                             gsl_matrix * q, gsl_matrix * r,
                             gsl_vector * tau,
                             gsl_permutation * p,
                             int *signum,
                             gsl_vector * norm);

GSL_FUN int gsl_linalg_QRPT_solve (const gsl_matrix * QR,
                           const gsl_vector * tau,
                           const gsl_permutation * p,
                           const gsl_vector * b,
                           gsl_vector * x);

GSL_FUN int gsl_linalg_QRPT_lssolve (const gsl_matrix * QR,
                             const gsl_vector * tau,
                             const gsl_permutation * p,
                             const gsl_vector * b,
                             gsl_vector * x,
                             gsl_vector * residual);

GSL_FUN int gsl_linalg_QRPT_lssolve2 (const gsl_matrix * QR,
                              const gsl_vector * tau,
                              const gsl_permutation * p,
                              const gsl_vector * b,
                              const size_t rank,
                              gsl_vector * x,
                              gsl_vector * residual);

GSL_FUN int gsl_linalg_QRPT_svx (const gsl_matrix * QR,
                         const gsl_vector * tau,
                         const gsl_permutation * p,
                         gsl_vector * x);

GSL_FUN int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,
                             const gsl_matrix * R,
                             const gsl_permutation * p,
                             const gsl_vector * b,
                             gsl_vector * x);

GSL_FUN int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
                             const gsl_permutation * p,
                             const gsl_vector * b,
                             gsl_vector * x);

GSL_FUN int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
                           const gsl_permutation * p,
                           gsl_vector * x);

GSL_FUN int gsl_linalg_QRPT_update (gsl_matrix * Q,
                            gsl_matrix * R,
                            const gsl_permutation * p,
                            gsl_vector * u,
                            const gsl_vector * v);

GSL_FUN size_t gsl_linalg_QRPT_rank (const gsl_matrix * QR, const double tol);

GSL_FUN int gsl_linalg_QRPT_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work);

/* triangle on top of rectangle QR decomposition */

GSL_FUN int gsl_linalg_QR_TR_decomp (gsl_matrix * S, gsl_matrix * A, gsl_matrix * T);

/* QL decomposition */

GSL_FUN int gsl_linalg_QL_decomp (gsl_matrix * A, gsl_vector * tau);

GSL_FUN int gsl_linalg_QL_unpack (const gsl_matrix * QL, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L);

/* COD decomposition */

GSL_FUN int gsl_linalg_COD_decomp(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
                          gsl_permutation * p, size_t * rank, gsl_vector * work);

GSL_FUN int gsl_linalg_COD_decomp_e(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
                            gsl_permutation * p, double tol, size_t * rank, gsl_vector * work);

GSL_FUN int gsl_linalg_COD_lssolve (const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
                            const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
                            gsl_vector * x, gsl_vector * residual);

GSL_FUN int
gsl_linalg_COD_lssolve2 (const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
                         const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
                         gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work);

GSL_FUN int gsl_linalg_COD_unpack(const gsl_matrix * QRZT, const gsl_vector * tau_Q,
                          const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q,
                          gsl_matrix * R, gsl_matrix * Z);

GSL_FUN int gsl_linalg_COD_matZ(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
                        gsl_matrix * A, gsl_vector * work);

/* LQ decomposition */

GSL_FUN int gsl_linalg_LQ_decomp (gsl_matrix * A, gsl_vector * tau);

GSL_FUN int gsl_linalg_LQ_lssolve(const gsl_matrix * LQ, const gsl_vector * tau,
                          const gsl_vector * b, gsl_vector * x, gsl_vector * residual);

GSL_FUN int gsl_linalg_LQ_QTvec(const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v);

GSL_FUN int gsl_linalg_LQ_solve_T (const gsl_matrix * LQ, const gsl_vector * tau,
			 const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_LQ_svx_T (const gsl_matrix * LQ, const gsl_vector * tau,
                         gsl_vector * x);

GSL_FUN int gsl_linalg_LQ_lssolve_T (const gsl_matrix * LQ, const gsl_vector * tau,
			   const gsl_vector * b, gsl_vector * x,
			   gsl_vector * residual);

GSL_FUN int gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b,
			  gsl_vector * x);

GSL_FUN int gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x);

GSL_FUN int gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b,
			gsl_vector * x);

GSL_FUN int gsl_linalg_LQ_vecQ (const gsl_matrix * LQ, const gsl_vector * tau,
			gsl_vector * v);

GSL_FUN int gsl_linalg_LQ_vecQT (const gsl_matrix * LQ, const gsl_vector * tau,
			 gsl_vector * v);

GSL_FUN int gsl_linalg_LQ_unpack (const gsl_matrix * LQ, const gsl_vector * tau,
			  gsl_matrix * Q, gsl_matrix * L);

GSL_FUN int gsl_linalg_LQ_update (gsl_matrix * Q, gsl_matrix * R,
			  const gsl_vector * v, gsl_vector * w);
GSL_FUN int gsl_linalg_LQ_LQsolve (gsl_matrix * Q, gsl_matrix * L,
			   const gsl_vector * b, gsl_vector * x);

/* P^T L Q decomposition */

GSL_FUN int gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau,
			    gsl_permutation * p, int *signum,
			    gsl_vector * norm);

GSL_FUN int gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q,
			     gsl_matrix * r, gsl_vector * tau,
			     gsl_permutation * p, int *signum,
			     gsl_vector * norm);

GSL_FUN int gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
			   const gsl_vector * tau,
			   const gsl_permutation * p,
			   const gsl_vector * b,
			   gsl_vector * x);

GSL_FUN int gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
                           const gsl_vector * tau,
                           const gsl_permutation * p,
                           gsl_vector * x);

GSL_FUN int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
			     const gsl_permutation * p,
			     const gsl_vector * b,
			     gsl_vector * x);

GSL_FUN int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
			    const gsl_permutation * p,
			    const gsl_vector * b,
			    gsl_vector * x);

GSL_FUN int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
			  const gsl_permutation * p,
			  gsl_vector * x);

GSL_FUN int gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
			    const gsl_permutation * p,
			    const gsl_vector * v, gsl_vector * w);

/* Cholesky Decomposition */

GSL_FUN int gsl_linalg_cholesky_decomp (gsl_matrix * A);
GSL_FUN int gsl_linalg_cholesky_decomp1 (gsl_matrix * A);

GSL_FUN int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky,
                               const gsl_vector * b,
                               gsl_vector * x);
GSL_FUN int gsl_linalg_cholesky_solve_mat (const gsl_matrix * cholesky,
                                   const gsl_matrix * B,
                                   gsl_matrix * X);

GSL_FUN int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky,
                             gsl_vector * x);
GSL_FUN int gsl_linalg_cholesky_svx_mat (const gsl_matrix * cholesky,
                                 gsl_matrix * X);

GSL_FUN int gsl_linalg_cholesky_invert(gsl_matrix * cholesky);

/* Cholesky decomposition with unit-diagonal triangular parts.
 *   A = L D L^T, where diag(L) = (1,1,...,1).
 *   Upon exit, A contains L and L^T as for Cholesky, and
 *   the diagonal of A is (1,1,...,1). The vector Dis set
 *   to the diagonal elements of the diagonal matrix D.
 */
GSL_FUN int gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D);

GSL_FUN int gsl_linalg_cholesky_scale(const gsl_matrix * A, gsl_vector * S);

GSL_FUN int gsl_linalg_cholesky_scale_apply(gsl_matrix * A, const gsl_vector * S);

GSL_FUN int gsl_linalg_cholesky_decomp2(gsl_matrix * A, gsl_vector * S);

GSL_FUN int gsl_linalg_cholesky_svx2 (const gsl_matrix * LLT,
                              const gsl_vector * S,
                              gsl_vector * x);

GSL_FUN int gsl_linalg_cholesky_solve2 (const gsl_matrix * LLT,
                                const gsl_vector * S,
                                const gsl_vector * b,
                                gsl_vector * x);

GSL_FUN int gsl_linalg_cholesky_rcond (const gsl_matrix * LLT, double * rcond,
                               gsl_vector * work);

/* Complex Cholesky Decomposition */

GSL_FUN int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A);

GSL_FUN int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                       const gsl_vector_complex * b,
                                       gsl_vector_complex * x);

GSL_FUN int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                     gsl_vector_complex * x);

GSL_FUN int gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * cholesky);

/* Pivoted Cholesky LDLT decomposition */

GSL_FUN int gsl_linalg_pcholesky_decomp (gsl_matrix * A, gsl_permutation * p);

GSL_FUN int gsl_linalg_pcholesky_solve(const gsl_matrix * LDLT,
                               const gsl_permutation * p,
                               const gsl_vector * b,
                               gsl_vector * x);

GSL_FUN int gsl_linalg_pcholesky_svx(const gsl_matrix * LDLT,
                             const gsl_permutation * p,
                             gsl_vector * x);

GSL_FUN int gsl_linalg_pcholesky_decomp2(gsl_matrix * A, gsl_permutation * p,
                                 gsl_vector * S);

GSL_FUN int gsl_linalg_pcholesky_solve2(const gsl_matrix * LDLT,
                                const gsl_permutation * p,
                                const gsl_vector * S,
                                const gsl_vector * b,
                                gsl_vector * x);

GSL_FUN int gsl_linalg_pcholesky_svx2(const gsl_matrix * LDLT,
                              const gsl_permutation * p,
                              const gsl_vector * S,
                              gsl_vector * x);

GSL_FUN int gsl_linalg_pcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
                                gsl_matrix * Ainv);

GSL_FUN int gsl_linalg_pcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
                                double * rcond, gsl_vector * work);

/* Modified Cholesky decomposition */

GSL_FUN int gsl_linalg_mcholesky_decomp (gsl_matrix * A, gsl_permutation * p,
                                 gsl_vector * E);

GSL_FUN int gsl_linalg_mcholesky_solve(const gsl_matrix * LDLT,
                               const gsl_permutation * p,
                               const gsl_vector * b,
                               gsl_vector * x);

GSL_FUN int gsl_linalg_mcholesky_svx(const gsl_matrix * LDLT,
                             const gsl_permutation * p,
                             gsl_vector * x);

GSL_FUN int gsl_linalg_mcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
                                double * rcond, gsl_vector * work);

GSL_FUN int gsl_linalg_mcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
                                gsl_matrix * Ainv);

/* Banded Cholesky decomposition */

GSL_FUN int gsl_linalg_cholesky_band_decomp(gsl_matrix * A);

GSL_FUN int gsl_linalg_cholesky_band_solve (const gsl_matrix * LLT, const gsl_vector * b,
                                    gsl_vector * x);

GSL_FUN int gsl_linalg_cholesky_band_svx (const gsl_matrix * LLT, gsl_vector * x);

GSL_FUN int gsl_linalg_cholesky_band_invert (const gsl_matrix * LLT, gsl_matrix * Ainv);

GSL_FUN int gsl_linalg_cholesky_band_unpack (const gsl_matrix * LLT, gsl_matrix * L);

GSL_FUN int gsl_linalg_cholesky_band_rcond (const gsl_matrix * LLT, double * rcond, gsl_vector * work);

/* L D L^T decomposition */

GSL_FUN int gsl_linalg_ldlt_decomp (gsl_matrix * A);

GSL_FUN int gsl_linalg_ldlt_solve (const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_ldlt_svx (const gsl_matrix * LDLT, gsl_vector * x);

GSL_FUN int gsl_linalg_ldlt_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work);

/* Banded L D L^T decomposition */

GSL_FUN int gsl_linalg_ldlt_band_decomp (gsl_matrix * A);

GSL_FUN int gsl_linalg_ldlt_band_solve (const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x);

GSL_FUN int gsl_linalg_ldlt_band_svx (const gsl_matrix * LDLT, gsl_vector * x);

GSL_FUN int gsl_linalg_ldlt_band_unpack (const gsl_matrix * LDLT, gsl_matrix * L, gsl_vector * D);

GSL_FUN int gsl_linalg_ldlt_band_rcond (const gsl_matrix * LDLT, double * rcond, gsl_vector * work);

/* Symmetric to symmetric tridiagonal decomposition */

GSL_FUN int gsl_linalg_symmtd_decomp (gsl_matrix * A,
                              gsl_vector * tau);

GSL_FUN int gsl_linalg_symmtd_unpack (const gsl_matrix * A,
                              const gsl_vector * tau,
                              gsl_matrix * Q,
                              gsl_vector * diag,
                              gsl_vector * subdiag);

GSL_FUN int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
                                gsl_vector * diag,
                                gsl_vector * subdiag);

/* Hermitian to symmetric tridiagonal decomposition */

GSL_FUN int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A,
                              gsl_vector_complex * tau);

GSL_FUN int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
                              const gsl_vector_complex * tau,
                              gsl_matrix_complex * U,
                              gsl_vector * diag,
                              gsl_vector * sudiag);

GSL_FUN int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
                                gsl_vector * diag,
                                gsl_vector * subdiag);

/* Linear Solve Using Householder Transformations

 * exceptions:
 */

GSL_FUN int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x);
GSL_FUN int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x);

/* Linear solve for a symmetric tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *     diag[0]  offdiag[0]             0    ...
 *  offdiag[0]     diag[1]    offdiag[1]    ...
 *           0  offdiag[1]       diag[2]    ...
 *           0           0    offdiag[2]    ...
 *         ...         ...           ...    ...
 */
GSL_FUN int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag,
                                   const gsl_vector * offdiag,
                                   const gsl_vector * b,
                                   gsl_vector * x);

/* Linear solve for a nonsymmetric tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *       diag[0]  abovediag[0]              0    ...
 *  belowdiag[0]       diag[1]   abovediag[1]    ...
 *             0  belowdiag[1]        diag[2]    ...
 *             0             0   belowdiag[2]    ...
 *           ...           ...            ...    ...
 */
GSL_FUN int gsl_linalg_solve_tridiag (const gsl_vector * diag,
                                   const gsl_vector * abovediag,
                                   const gsl_vector * belowdiag,
                                   const gsl_vector * b,
                                   gsl_vector * x);

/* Linear solve for a symmetric cyclic tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
 *   offdiag[0]     diag[1]    offdiag[1]   .....
 *            0  offdiag[1]       diag[2]   .....
 *            0           0    offdiag[2]   .....
 *          ...         ...
 * offdiag[N-1]         ...
 */
GSL_FUN int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag,
                                       const gsl_vector * offdiag,
                                       const gsl_vector * b,
                                       gsl_vector * x);

/* Linear solve for a nonsymmetric cyclic tridiagonal system.

 * The input vectors represent the NxN matrix as follows:
 *
 *        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
 *   belowdiag[0]       diag[1]  abovediag[1]   .....
 *              0  belowdiag[1]       diag[2]
 *              0             0  belowdiag[2]   .....
 *            ...           ...
 * abovediag[N-1]           ...
 */
GSL_FUN int gsl_linalg_solve_cyc_tridiag (const gsl_vector * diag,
                                  const gsl_vector * abovediag,
                                  const gsl_vector * belowdiag,
                                  const gsl_vector * b,
                                  gsl_vector * x);

/* Bidiagonal decomposition */

GSL_FUN int gsl_linalg_bidiag_decomp (gsl_matrix * A,
                              gsl_vector * tau_U,
                              gsl_vector * tau_V);

GSL_FUN int gsl_linalg_bidiag_unpack (const gsl_matrix * A,
                              const gsl_vector * tau_U,
                              gsl_matrix * U,
                              const gsl_vector * tau_V,
                              gsl_matrix * V,
                              gsl_vector * diag,
                              gsl_vector * superdiag);

GSL_FUN int gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
                               gsl_vector * tau_U,
                               gsl_vector * tau_V,
                               gsl_matrix * V);

GSL_FUN int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
                                gsl_vector * diag,
                                gsl_vector * superdiag);

/* Balancing */

GSL_FUN int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D);
GSL_FUN int gsl_linalg_balance_accum (gsl_matrix * A, gsl_vector * D);
GSL_FUN int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D);

/* condition estimation */

GSL_FUN int gsl_linalg_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work);
GSL_FUN int gsl_linalg_tri_upper_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work);
GSL_FUN int gsl_linalg_tri_lower_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work);
GSL_FUN int gsl_linalg_invnorm1(const size_t N,
                        int (* Ainvx)(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params),
                        void * params, double * Ainvnorm, gsl_vector * work);

/* triangular matrices */

GSL_FUN int gsl_linalg_tri_upper_invert(gsl_matrix * T);
GSL_FUN int gsl_linalg_tri_lower_invert(gsl_matrix * T);
GSL_FUN int gsl_linalg_tri_upper_unit_invert(gsl_matrix * T);
GSL_FUN int gsl_linalg_tri_lower_unit_invert(gsl_matrix * T);

GSL_FUN int gsl_linalg_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T);
GSL_FUN int gsl_linalg_complex_tri_invert(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T);

GSL_FUN int gsl_linalg_tri_LTL(gsl_matrix * L);
GSL_FUN int gsl_linalg_tri_UL(gsl_matrix * LU);
GSL_FUN int gsl_linalg_complex_tri_LHL(gsl_matrix_complex * L);
GSL_FUN int gsl_linalg_complex_tri_UL(gsl_matrix_complex * LU);

GSL_FUN INLINE_DECL void gsl_linalg_givens (const double a, const double b,
                                    double *c, double *s);
GSL_FUN INLINE_DECL void gsl_linalg_givens_gv (gsl_vector * v, const size_t i,
                                       const size_t j, const double c,
                                       const double s);

#ifdef HAVE_INLINE

/* Generate a Givens rotation (cos,sin) which takes v=(x,y) to (|v|,0)
   From Golub and Van Loan, "Matrix Computations", Section 5.1.8 */
INLINE_FUN
void
gsl_linalg_givens (const double a, const double b, double *c, double *s)
{
  if (b == 0)
    {
      *c = 1;
      *s = 0;
    }
  else if (fabs (b) > fabs (a))
    {
      double t = -a / b;
      double s1 = 1.0 / sqrt (1 + t * t);
      *s = s1;
      *c = s1 * t;
    }
  else
    {
      double t = -b / a;
      double c1 = 1.0 / sqrt (1 + t * t);
      *c = c1;
      *s = c1 * t;
    }
} /* gsl_linalg_givens() */

INLINE_FUN
void
gsl_linalg_givens_gv (gsl_vector * v, const size_t i, const size_t j,
                      const double c, const double s)
{
  /* Apply rotation to vector v' = G^T v */

  double vi = gsl_vector_get (v, i);
  double vj = gsl_vector_get (v, j);
  gsl_vector_set (v, i, c * vi - s * vj);
  gsl_vector_set (v, j, s * vi + c * vj);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_LINALG_H__ */


text:
/*
 * Copyright 2017 RIS.
 * Authored by Ron Dahlgren <ron@red83.net>
 */
#include "entities.h"
#include "utility/lists.h"
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

vector* entity_vec;

entity_info entity_template;

entity_info** working_copy;

static void no_op_renderer(double deltaT, entity_info *e) { }

static sd_result initialize() {
  gsl_matrix_view frame_view;
  const int starting_size = 10;
  working_copy = malloc(sizeof(void*) * starting_size);

  entity_vec = vector_ops.vector_alloc(starting_size);

  /* Set the desired template values */
  entity_template.shader = shaders.load_material(MAT_ID_BASIC_ALPHA);
  entity_template.renderer = no_op_renderer;
  bzero(&entity_template.animation_state, sizeof(animation_state));
  bzero(&entity_template.math_context, sizeof(math_context));
  /* First the reference frame */
  frame_view = gsl_matrix_view_array(entity_template.reference_frame, 4, 4);
  gsl_matrix_set_identity(&frame_view.matrix);
  /* Now the transform */
  frame_view = gsl_matrix_view_array(entity_template.transform, 4, 4);
  gsl_matrix_set_identity(&frame_view.matrix);

  entity_template.initialize = NULL;
  entity_template.terminate = NULL;
  entity_template.pre_update = NULL;
  entity_template.post_update = NULL;

  entity_template.data = NULL;

  return SD_OK;
}

static sd_result release_resources() {
  entity_info *cursor = NULL;
  while ((cursor = vector_ops.pop(entity_vec)) != NULL) {
    if (cursor->terminate != NULL) { cursor->terminate(cursor); }
    entities.entity_free(cursor);
  }
  return SD_OK;
}

static sd_result terminate() {
  release_resources();
  vector_ops.vector_free(&entity_vec);
  free(working_copy);

  return SD_OK;
}

void add_entity(entity_info* info) {
  vector_ops.append(entity_vec, info);
  // Resize buffer
  int size = entity_vec->size;
  free(working_copy);
  working_copy = malloc(sizeof(void*) * size);
  // initialize the entity
  if (info->initialize != NULL) { info->initialize(info); }
}

void remove_entity(entity_info* info) {
  int located_index = -1;
  unsigned int i;

  for (i = 0; i < entity_vec->size; ++i) {
    if (entity_vec->storage[i] == info) {
      located_index = i;
      break;
    }
  }
  if (located_index >= 0) {
    vector_ops.remove(entity_vec, located_index);
    // terminate the entity
    if (info->terminate != NULL) { info->terminate(info); }
  }
}

void update(double deltaT) {
  int size = entity_vec->size;
  /* TODO Safety lock here */
  memcpy(working_copy, entity_vec->storage, sizeof(void*) * size);
  /* TODO End safety lock here */

  /* Update the entities */
  entity_info *cursor = NULL;
  for (int i = 0; i < size; i++) {
    cursor = working_copy[i];
    if (cursor->pre_update) { cursor->pre_update(deltaT, cursor); }
    // UPDATE!
    // animation
    transform_animation.animate_entity(deltaT, working_copy[i]);

    gsl_matrix_view transform_view = gsl_matrix_view_array(
        cursor->transform, 4, 4);

    gsl_matrix_view frame_view = gsl_matrix_view_array(
        cursor->reference_frame, 4, 4);

    // Apply the animation if necessary
    if (cursor->animation_state.current_cycle != NULL) {
      animation_state *state = &cursor->animation_state;
      int reverse = state->playback_speed < 0;

      state->current_cycle->get_transform(
          state->current_cycle_time,
          state->playback_mode,
          reverse,
          state->current_cycle->data,
          state->animation_solution);

      gsl_matrix_view animation_view = gsl_matrix_view_array(
          state->animation_solution, 4, 4);

      // T = Frame * Anim
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &frame_view.matrix, &animation_view.matrix, 0.0, &transform_view.matrix);
    } else {
      gsl_matrix_memcpy(&transform_view.matrix, &frame_view.matrix);
    }
    // physics?

    if (cursor->post_update) { cursor->post_update(deltaT, cursor); }
  }
}

entity_info* entity_alloc() {
  entity_info* result = (entity_info*)malloc(sizeof(entity_info));
  memcpy(result, &entity_template, sizeof(entity_info));
  return result;
}
void entity_free(entity_info* info) { free(info); }

entities_struct const entities = {
  .initialize=initialize,
  .release_resources=release_resources,
  .terminate=terminate,
  .entity_alloc=entity_alloc,
  .entity_free=entity_free,
  .add_entity=add_entity,
  .remove_entity=remove_entity,
  .update=update
};

ext: c


text:
static char help[] =
"Solves the p-Helmholtz equation in 2D using Q_1 FEM.  Option prefix -ph_.\n"
"Problem is posed as minimizing this objective functional over W^{1,p}\n"
"for p>1:\n"
"    I[u] = int_Omega (1/p) |grad u|^p + (1/2) u^2 - f u.\n"
"The strong form equation, namely setting the gradient to zero, is a PDE\n"
"    - div( |grad u|^{p-2} grad u ) + u = f\n"
"subject to homogeneous Neumann boundary conditions.  Implements objective\n"
"and gradient (residual) but no Hessian (Jacobian).  Defaults to linear\n"
"problem (p=2) and quadrature degree 2.  Can be run with only an objective\n" "function; use -ph_no_gradient -snes_fd_function.\n\n";

#include <petsc.h>
#include "../interlude/quadrature.h"

typedef struct {
    PetscReal  p, eps;
    PetscInt   quadpts;
    PetscReal  (*f)(PetscReal x, PetscReal y, PetscReal p, PetscReal eps);
} PHelmCtx;

static PetscReal f_constant(PetscReal x, PetscReal y, PetscReal p, PetscReal eps) {
    return 1.0;
}

static PetscReal u_exact_cosines(PetscReal x, PetscReal y, PetscReal p, PetscReal eps) {
    return PetscCosReal(PETSC_PI * x) * PetscCosReal(PETSC_PI * y);
}

static PetscReal f_cosines(PetscReal x, PetscReal y, PetscReal p, PetscReal eps) {
    const PetscReal uu = u_exact_cosines(x,y,p,eps),
                    pi2 = PETSC_PI * PETSC_PI,
                    lapu = - 2 * pi2 * uu;
    if (p == 2.0) {
        return - lapu + uu;
    } else {
        const PetscReal
            ux = - PETSC_PI * PetscSinReal(PETSC_PI * x)
                 * PetscCosReal(PETSC_PI * y),
            uy = - PETSC_PI * PetscCosReal(PETSC_PI * x)
                 * PetscSinReal(PETSC_PI * y),
            // note regularization changes f(x,y) but not u(x,y):
            w = ux * ux + uy * uy + eps * eps,
            pi3 = pi2 * PETSC_PI,
            wx = pi3 * PetscSinReal(2 * PETSC_PI * x)
                 * PetscCosReal(2 * PETSC_PI * y),
            wy = pi3 * PetscCosReal(2 * PETSC_PI * x)
                 * PetscSinReal(2 * PETSC_PI * y);
        const PetscReal s = (p - 2) / 2;  //  -1/2 <= s <= 0
        return - s * PetscPowScalar(w,s-1) * (wx * ux + wy * uy)
               - PetscPowScalar(w,s) * lapu + uu;
    }
}

typedef enum {CONSTANT, COSINES} ProblemType;
static const char* ProblemTypes[] = {"constant","cosines",
                                     "ProblemType", "", NULL};

extern PetscErrorCode GetVecFromFunction(DMDALocalInfo*, Vec,
                         PetscReal (*)(PetscReal, PetscReal, PetscReal, PetscReal), PHelmCtx*);
extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*, PetscReal**, PetscReal*, PHelmCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, PetscReal**, PetscReal**, PHelmCtx*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da;
    SNES           snes;
    Vec            u_initial, u, u_exact;
    PHelmCtx       user;
    DMDALocalInfo  info;
    ProblemType    problem = COSINES;
    PetscBool      no_objective = PETSC_FALSE,
                   no_gradient = PETSC_FALSE,
                   exact_init = PETSC_FALSE,
                   view_f = PETSC_FALSE;
    PetscReal      err;

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

    user.p = 2.0;
    user.eps = 0.0;
    user.quadpts = 2;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"ph_",
                  "p-Helmholtz solver options",""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-eps",
                  "regularization parameter eps",
                  "phelm.c",user.eps,&(user.eps),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-exact_init",
                  "use exact solution to initialize",
                  "phelm.c",exact_init,&(exact_init),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-no_objective",
                  "do not set the objective evaluation function",
                  "phelm.c",no_objective,&(no_objective),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-no_gradient",
                  "do not set the residual evaluation function",
                  "phelm.c",no_gradient,&(no_gradient),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-p",
                  "exponent p > 1",
                  "phelm.c",user.p,&(user.p),NULL); CHKERRQ(ierr);
    if (user.p < 1.0) {
         SETERRQ(PETSC_COMM_SELF,1,"p >= 1 required");
    }
    if (user.p == 1.0) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,
            "WARNING: well-posedness only known for p > 1\n"); CHKERRQ(ierr);
    }
    ierr = PetscOptionsEnum("-problem",
                  "problem type determines right side f(x,y)",
                  "phelm.c",ProblemTypes,(PetscEnum)problem,(PetscEnum*)&problem,
                  NULL); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-quadpts",
                  "number n of quadrature points in each direction (= 1,2,3 only)",
                  "phelm.c",user.quadpts,&(user.quadpts),NULL); CHKERRQ(ierr);
    if ((user.quadpts < 1) || (user.quadpts > 3)) {
        SETERRQ(PETSC_COMM_SELF,3,"quadrature points n=1,2,3 only");
    }
    ierr = PetscOptionsBool("-view_f",
                  "view right-hand side to STDOUT",
                  "phelm.c",view_f,&(view_f),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    ierr = DMDACreate2d(PETSC_COMM_WORLD,
           DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
           2,2,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);
    ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,-1.0,-1.0); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
    if (!no_objective) {
        ierr = DMDASNESSetObjectiveLocal(da,
             (DMDASNESObjective)FormObjectiveLocal,&user); CHKERRQ(ierr);
    }
    if (no_gradient) {
        // why isn't this the default?  why no programmatic way to set?
        ierr = PetscOptionsSetValue(NULL,"-snes_fd_function_eps","0.0"); CHKERRQ(ierr);
    } else {
        ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
    }
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

    // set initial iterate and right-hand side
    ierr = DMCreateGlobalVector(da,&u_initial);CHKERRQ(ierr);
    ierr = VecSet(u_initial,0.5); CHKERRQ(ierr);
    switch (problem) {
        case CONSTANT:
            if (exact_init) {
                ierr = VecSet(u_initial,1.0); CHKERRQ(ierr);
            }
            user.f = &f_constant;
            break;
        case COSINES:
            if (exact_init) {
                ierr = GetVecFromFunction(&info,u_initial,&u_exact_cosines,&user); CHKERRQ(ierr);
            }
            user.f = &f_cosines;
            break;
        default:
            SETERRQ(PETSC_COMM_SELF,4,"unknown problem type\n");
    }

    // optionally view right-hand-side on initial grid
    if (view_f) {
        Vec vf;
        ierr = VecDuplicate(u_initial,&vf); CHKERRQ(ierr);
        switch (problem) {
            case CONSTANT:
                ierr = VecSet(vf,1.0); CHKERRQ(ierr);
                break;
            case COSINES:
                ierr = GetVecFromFunction(&info,vf,&f_cosines,&user); CHKERRQ(ierr);
                break;
        }
        ierr = VecView(vf,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
        VecDestroy(&vf);
    }

    // solve and clean up
    ierr = SNESSolve(snes,NULL,u_initial); CHKERRQ(ierr);
    ierr = VecDestroy(&u_initial); CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);
    ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);

    // evaluate numerical error
    ierr = VecDuplicate(u,&u_exact); CHKERRQ(ierr);
    switch (problem) {
        case CONSTANT:
            ierr = VecSet(u_exact,1.0); CHKERRQ(ierr);
            break;
        case COSINES:
            ierr = GetVecFromFunction(&info,u_exact,&u_exact_cosines,&user); CHKERRQ(ierr);
            break;
    }
    ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr);    // u <- u + (-1.0) uexact
    ierr = VecNorm(u,NORM_INFINITY,&err); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
        "done on %d x %d grid with p=%.3f ...\n"
        "  numerical error:  |u-u_exact|_inf = %.3e\n",
        info.mx,info.my,user.p,err); CHKERRQ(ierr);

    VecDestroy(&u_exact);  SNESDestroy(&snes);
    return PetscFinalize();
}

PetscErrorCode GetVecFromFunction(DMDALocalInfo *info, Vec w,
         PetscReal (*fcn)(PetscReal, PetscReal, PetscReal, PetscReal),
         PHelmCtx *user) {
    PetscErrorCode  ierr;
    const PetscReal hx = 1.0 / (info->mx - 1), hy = 1.0 / (info->my - 1);
    PetscReal       x, y, **aw;
    PetscInt        i, j;
    ierr = DMDAVecGetArray(info->da,w,&aw); CHKERRQ(ierr);
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = i * hx;
            aw[j][i] = (*fcn)(x,y,user->p,user->eps);
        }
    }
    ierr = DMDAVecRestoreArray(info->da,w,&aw); CHKERRQ(ierr);
    return 0;
}

//STARTFEM
static PetscReal xiL[4]  = { 1.0, -1.0, -1.0,  1.0},
                 etaL[4] = { 1.0,  1.0, -1.0, -1.0};

static PetscReal chi(PetscInt L, PetscReal xi, PetscReal eta) {
    return 0.25 * (1.0 + xiL[L] * xi) * (1.0 + etaL[L] * eta);
}

// evaluate v(xi,eta) on reference element using local node numbering
static PetscReal eval(const PetscReal v[4], PetscReal xi, PetscReal eta) {
    return   v[0] * chi(0,xi,eta) + v[1] * chi(1,xi,eta)
           + v[2] * chi(2,xi,eta) + v[3] * chi(3,xi,eta);
}

typedef struct {
    PetscReal  xi, eta;
} gradRef;

static gradRef dchi(PetscInt L, PetscReal xi, PetscReal eta) {
    const gradRef result = {0.25 * xiL[L]  * (1.0 + etaL[L] * eta),
                            0.25 * etaL[L] * (1.0 + xiL[L]  * xi)};
    return result;
}

// evaluate partial derivs of v(xi,eta) on reference element
static gradRef deval(const PetscReal v[4], PetscReal xi, PetscReal eta) {
    gradRef   sum = {0.0,0.0}, tmp;
    PetscInt  L;
    for (L=0; L<4; L++) {
        tmp = dchi(L,xi,eta);
        sum.xi += v[L] * tmp.xi;  sum.eta += v[L] * tmp.eta;
    }
    return sum;
}

static PetscReal GradInnerProd(PetscReal hx, PetscReal hy,
                               gradRef du, gradRef dv) {
    const PetscReal cx = 4.0 / (hx * hx),  cy = 4.0 / (hy * hy);
    return cx * du.xi * dv.xi + cy * du.eta * dv.eta;
}

static PetscReal GradPow(PetscReal hx, PetscReal hy,
                         gradRef du, PetscReal P, PetscReal eps) {
    return PetscPowScalar(GradInnerProd(hx,hy,du,du) + eps*eps, P/2.0);
}
//ENDFEM

/* FLOPS:  (counting PetscPowScalar as 1)
     chi = 6
     eval = 4*6+7 = 31
     dchi = 8
     deval = 4*8+4 = 36
     GradInnerProd = 9
     GradPow = 9+4 = 13
     ObjIntegrandRef = deval + 2*eval + GradPow + 10 = 121
     FunIntegrandRef = chi + dchi + 2*eval + deval + GradPo + GradInnerProd + 9
                     = 143
*/

//STARTOBJECTIVE
static PetscReal ObjIntegrandRef(DMDALocalInfo *info,
                       const PetscReal ff[4], const PetscReal uu[4],
                       PetscReal xi, PetscReal eta, PHelmCtx *user) {
    const gradRef    du = deval(uu,xi,eta);
    const PetscReal  hx = 1.0 / (info->mx-1),  hy = 1.0 / (info->my-1),
                     u = eval(uu,xi,eta);
    return GradPow(hx,hy,du,user->p,0.0) / user->p + 0.5 * u * u
           - eval(ff,xi,eta) * u;
}

PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info, PetscReal **au,
                                  PetscReal *obj, PHelmCtx *user) {
  PetscErrorCode  ierr;
  const PetscReal hx = 1.0 / (info->mx-1),  hy = 1.0 / (info->my-1);
  const Quad1D    q = gausslegendre[user->quadpts-1];
  PetscReal       x, y, lobj = 0.0;
  PetscInt        i,j,r,s;
  MPI_Comm        com;

  // loop over all elements
  for (j = info->ys; j < info->ys + info->ym; j++) {
      if (j == 0)
          continue;
      y = j * hy;
      for (i = info->xs; i < info->xs + info->xm; i++) {
          if (i == 0)
              continue;
          x = i * hx;
          const PetscReal ff[4] = {user->f(x,y,user->p,user->eps),
                                   user->f(x-hx,y,user->p,user->eps),
                                   user->f(x-hx,y-hy,user->p,user->eps),
                                   user->f(x,y-hy,user->p,user->eps)};
          const PetscReal uu[4] = {au[j][i],au[j][i-1],
                                   au[j-1][i-1],au[j-1][i]};
          // loop over quadrature points on this element
          for (r = 0; r < q.n; r++) {
              for (s = 0; s < q.n; s++) {
                  lobj += q.w[r] * q.w[s]
                          * ObjIntegrandRef(info,ff,uu,
                                            q.xi[r],q.xi[s],user);
              }
          }
      }
  }
  lobj *= hx * hy / 4.0;  // from change of variables formula
  ierr = PetscObjectGetComm((PetscObject)(info->da),&com); CHKERRQ(ierr);
  ierr = MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,com); CHKERRQ(ierr);
  ierr = PetscLogFlops(129*info->xm*info->ym); CHKERRQ(ierr);
  return 0;
}
//ENDOBJECTIVE

//STARTFUNCTION
static PetscReal IntegrandRef(DMDALocalInfo *info, PetscInt L,
                     const PetscReal ff[4], const PetscReal uu[4],
                     PetscReal xi, PetscReal eta, PHelmCtx *user) {
  const gradRef    du    = deval(uu,xi,eta),
                   dchiL = dchi(L,xi,eta);
  const PetscReal  hx = 1.0 / (info->mx-1),  hy = 1.0 / (info->my-1);
  return GradPow(hx,hy,du,user->p - 2.0,user->eps)
           * GradInnerProd(hx,hy,du,dchiL)
         + (eval(uu,xi,eta) - eval(ff,xi,eta)) * chi(L,xi,eta);
}

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscReal **au,
                                 PetscReal **FF, PHelmCtx *user) {
  PetscErrorCode ierr;
  const PetscReal hx = 1.0 / (info->mx-1),  hy = 1.0 / (info->my-1);
  const Quad1D    q = gausslegendre[user->quadpts-1];
  const PetscInt  li[4] = {0,-1,-1,0},  lj[4] = {0,0,-1,-1};
  PetscReal       x, y;
  PetscInt        i,j,l,r,s,PP,QQ;

  // clear residuals
  for (j = info->ys; j < info->ys + info->ym; j++)
      for (i = info->xs; i < info->xs + info->xm; i++)
          FF[j][i] = 0.0;

  // loop over all elements
  for (j = info->ys; j <= info->ys + info->ym; j++) {
      if ((j == 0) || (j > info->my-1))
          continue;
      y = j * hy;
      for (i = info->xs; i <= info->xs + info->xm; i++) {
          if ((i == 0) || (i > info->mx-1))
              continue;
          x = i * hx;
          const PetscReal ff[4] = {user->f(x,y,user->p,user->eps),
                                   user->f(x-hx,y,user->p,user->eps),
                                   user->f(x-hx,y-hy,user->p,user->eps),
                                   user->f(x,y-hy,user->p,user->eps)};
          const PetscReal uu[4] = {au[j][i],au[j][i-1],
                                   au[j-1][i-1],au[j-1][i]};
          // loop over corners of element i,j
          for (l = 0; l < 4; l++) {
              PP = i + li[l];
              QQ = j + lj[l];
              // only update residual if we own node
              if (PP >= info->xs && PP < info->xs + info->xm
                  && QQ >= info->ys && QQ < info->ys + info->ym) {
                  // loop over quadrature points
                  for (r = 0; r < q.n; r++) {
                      for (s = 0; s < q.n; s++) {
                         FF[QQ][PP]
                             += 0.25 * hx * hy * q.w[r] * q.w[s]
                                * IntegrandRef(info,l,ff,uu,
                                               q.xi[r],q.xi[s],user);
                      }
                  }
              }
          }
      }
  }
  ierr = PetscLogFlops((5+q.n*q.n*149)*(info->xm+1)*(info->ym+1)); CHKERRQ(ierr);
  return 0;
}
//ENDFUNCTION


text:
/* Author: G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>

/* Compile all the inline functions */

#define COMPILE_INLINE_STATIC
#include "build.h"
#include <gsl/gsl_qrng.h>

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <ceed.h>
#include <petsc.h>
#include "../include/bcfunctions.h"

// -----------------------------------------------------------------------------
// Boundary Conditions
// -----------------------------------------------------------------------------

// Diff boundary condition function
PetscErrorCode BCsDiff(PetscInt dim, PetscReal time, const PetscReal x[],
                       PetscInt num_comp_u, PetscScalar *u, void *ctx) {
  // *INDENT-OFF*
  #ifndef M_PI
  #define M_PI    3.14159265358979323846
  #endif
  // *INDENT-ON*
  const CeedScalar c[3] = { 0, 1., 2. };
  const CeedScalar k[3] = { 1., 2., 3. };

  PetscFunctionBeginUser;

  for (PetscInt i = 0; i < num_comp_u; i++)
    u[i] = sin(M_PI*(c[0] + k[0]*x[0])) *
           sin(M_PI*(c[1] + k[1]*x[1])) *
           sin(M_PI*(c[2] + k[2]*x[2]));

  PetscFunctionReturn(0);
}

// Mass boundary condition function
PetscErrorCode BCsMass(PetscInt dim, PetscReal time, const PetscReal x[],
                       PetscInt num_comp_u, PetscScalar *u, void *ctx) {
  PetscFunctionBeginUser;

  for (PetscInt i = 0; i < num_comp_u; i++)
    u[i] = PetscSqrtScalar(PetscSqr(x[0]) + PetscSqr(x[1]) +
                           PetscSqr(x[2]));

  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
/* matrix/gsl_matrix_complex_float.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_COMPLEX_FLOAT_H__
#define __GSL_MATRIX_COMPLEX_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_float.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  float * data;
  gsl_block_complex_float * block;
  int owner;
} gsl_matrix_complex_float ;

typedef struct
{
  gsl_matrix_complex_float matrix;
} _gsl_matrix_complex_float_view;

typedef _gsl_matrix_complex_float_view gsl_matrix_complex_float_view;

typedef struct
{
  gsl_matrix_complex_float matrix;
} _gsl_matrix_complex_float_const_view;

typedef const _gsl_matrix_complex_float_const_view gsl_matrix_complex_float_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_complex_float *
gsl_matrix_complex_float_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc_from_block (gsl_block_complex_float * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

GSL_EXPORT
gsl_matrix_complex_float *
gsl_matrix_complex_float_alloc_from_matrix (gsl_matrix_complex_float * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

GSL_EXPORT
gsl_vector_complex_float *
gsl_vector_complex_float_alloc_row_from_matrix (gsl_matrix_complex_float * m,
                                                const size_t i);

GSL_EXPORT
gsl_vector_complex_float *
gsl_vector_complex_float_alloc_col_from_matrix (gsl_matrix_complex_float * m,
                                                const size_t j);

GSL_EXPORT void gsl_matrix_complex_float_free (gsl_matrix_complex_float * m);

/* Views */

GSL_EXPORT
_gsl_matrix_complex_float_view
gsl_matrix_complex_float_submatrix (gsl_matrix_complex_float * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_complex_float_view
gsl_matrix_complex_float_row (gsl_matrix_complex_float * m, const size_t i);

GSL_EXPORT
_gsl_vector_complex_float_view
gsl_matrix_complex_float_column (gsl_matrix_complex_float * m, const size_t j);

GSL_EXPORT
_gsl_vector_complex_float_view
gsl_matrix_complex_float_diagonal (gsl_matrix_complex_float * m);

GSL_EXPORT
_gsl_vector_complex_float_view
gsl_matrix_complex_float_subdiagonal (gsl_matrix_complex_float * m, const size_t k);

GSL_EXPORT
_gsl_vector_complex_float_view
gsl_matrix_complex_float_superdiagonal (gsl_matrix_complex_float * m, const size_t k);

GSL_EXPORT
_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_array (float * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_array_with_tda (float * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_vector (gsl_vector_complex_float * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_float_view
gsl_matrix_complex_float_view_vector_with_tda (gsl_vector_complex_float * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_submatrix (const gsl_matrix_complex_float * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_row (const gsl_matrix_complex_float * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_column (const gsl_matrix_complex_float * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_diagonal (const gsl_matrix_complex_float * m);

GSL_EXPORT
_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_subdiagonal (const gsl_matrix_complex_float * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_complex_float_const_view
gsl_matrix_complex_float_const_superdiagonal (const gsl_matrix_complex_float * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_array (const float * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_array_with_tda (const float * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_vector (const gsl_vector_complex_float * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_complex_float_const_view
gsl_matrix_complex_float_const_view_vector_with_tda (const gsl_vector_complex_float * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT gsl_complex_float gsl_matrix_complex_float_get(const gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_EXPORT void gsl_matrix_complex_float_set(gsl_matrix_complex_float * m, const size_t i, const size_t j, const gsl_complex_float x);

GSL_EXPORT gsl_complex_float * gsl_matrix_complex_float_ptr(gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_EXPORT const gsl_complex_float * gsl_matrix_complex_float_const_ptr(const gsl_matrix_complex_float * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_complex_float_set_zero (gsl_matrix_complex_float * m);
GSL_EXPORT void gsl_matrix_complex_float_set_identity (gsl_matrix_complex_float * m);
GSL_EXPORT void gsl_matrix_complex_float_set_all (gsl_matrix_complex_float * m, gsl_complex_float x);

GSL_EXPORT int gsl_matrix_complex_float_fread (FILE * stream, gsl_matrix_complex_float * m) ;
GSL_EXPORT int gsl_matrix_complex_float_fwrite (FILE * stream, const gsl_matrix_complex_float * m) ;
GSL_EXPORT int gsl_matrix_complex_float_fscanf (FILE * stream, gsl_matrix_complex_float * m);
GSL_EXPORT int gsl_matrix_complex_float_fprintf (FILE * stream, const gsl_matrix_complex_float * m, const char * format);

GSL_EXPORT int gsl_matrix_complex_float_memcpy(gsl_matrix_complex_float * dest, const gsl_matrix_complex_float * src);
GSL_EXPORT int gsl_matrix_complex_float_swap(gsl_matrix_complex_float * m1, gsl_matrix_complex_float * m2);

GSL_EXPORT int gsl_matrix_complex_float_swap_rows(gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_complex_float_swap_columns(gsl_matrix_complex_float * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_complex_float_swap_rowcol(gsl_matrix_complex_float * m, const size_t i, const size_t j);

GSL_EXPORT int gsl_matrix_complex_float_transpose (gsl_matrix_complex_float * m);
GSL_EXPORT int gsl_matrix_complex_float_transpose_memcpy (gsl_matrix_complex_float * dest, const gsl_matrix_complex_float * src);

GSL_EXPORT int gsl_matrix_complex_float_isnull (const gsl_matrix_complex_float * m);

GSL_EXPORT int gsl_matrix_complex_float_add (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_EXPORT int gsl_matrix_complex_float_sub (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_EXPORT int gsl_matrix_complex_float_mul_elements (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_EXPORT int gsl_matrix_complex_float_div_elements (gsl_matrix_complex_float * a, const gsl_matrix_complex_float * b);
GSL_EXPORT int gsl_matrix_complex_float_scale (gsl_matrix_complex_float * a, const gsl_complex_float x);
GSL_EXPORT int gsl_matrix_complex_float_add_constant (gsl_matrix_complex_float * a, const gsl_complex_float x);
GSL_EXPORT int gsl_matrix_complex_float_add_diagonal (gsl_matrix_complex_float * a, const gsl_complex_float x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_complex_float_get_row(gsl_vector_complex_float * v, const gsl_matrix_complex_float * m, const size_t i);
GSL_EXPORT int gsl_matrix_complex_float_get_col(gsl_vector_complex_float * v, const gsl_matrix_complex_float * m, const size_t j);
GSL_EXPORT int gsl_matrix_complex_float_set_row(gsl_matrix_complex_float * m, const size_t i, const gsl_vector_complex_float * v);
GSL_EXPORT int gsl_matrix_complex_float_set_col(gsl_matrix_complex_float * m, const size_t j, const gsl_vector_complex_float * v);

#ifdef HAVE_INLINE

extern inline
gsl_complex_float
gsl_matrix_complex_float_get(const gsl_matrix_complex_float * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  gsl_complex_float zero = {{0,0}};

  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
    }
#endif
  return *(gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

extern inline
void
gsl_matrix_complex_float_set(gsl_matrix_complex_float * m,
                     const size_t i, const size_t j, const gsl_complex_float x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  *(gsl_complex_float *)(m->data + 2*(i * m->tda + j)) = x ;
}

extern inline
gsl_complex_float *
gsl_matrix_complex_float_ptr(gsl_matrix_complex_float * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

extern inline
const gsl_complex_float *
gsl_matrix_complex_float_const_ptr(const gsl_matrix_complex_float * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const gsl_complex_float *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_FLOAT_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#pragma once
#include "iarg.h"
#include "optioninfo.h"
#include "configaccess.h"
#include "format.h"
#include "streamreader.h"
#include <gsl/gsl>
#include <cmdlime/errors.h>
#include <cmdlime/customnames.h>
#include <sstream>
#include <functional>
#include <memory>

namespace cmdlime::detail{

template <typename T>
class Arg : public IArg{
public:
    Arg(std::string name,
        std::string type,
        std::function<T&()> argGetter)
        : info_(std::move(name), {}, std::move(type))
        , argGetter_(std::move(argGetter))
    {
    }

    OptionInfo& info() override
    {
        return info_;
    }

    const OptionInfo& info() const override
    {
        return info_;
    }

private:
    bool read(const std::string& data) override
    {
        auto stream = std::stringstream{data};
        return readFromStream(stream, argGetter_());
    }

private:
    OptionInfo info_;
    std::function<T&()> argGetter_;
};

template <>
inline bool Arg<std::string>::read(const std::string& data)
{
    argGetter_() = data;
    return true;
}

template<typename T, typename TConfig>
class ArgCreator{
    using NameProvider = typename Format<ConfigAccess<TConfig>::format()>::nameProvider;
public:
    ArgCreator(TConfig& cfg,
               const std::string& varName,
               const std::string& type,
               std::function<T&()> argGetter)
        : cfg_(cfg)
    {
        Expects(!varName.empty());
        Expects(!type.empty());
        arg_ = std::make_unique<Arg<T>>(NameProvider::fullName(varName),
                                        NameProvider::valueName(type),
                                        std::move(argGetter));
    }

    ArgCreator<T, TConfig>& operator<<(const std::string& info)
    {
        arg_->info().addDescription(info);
        return *this;
    }

    ArgCreator<T, TConfig>& operator<<(const Name& customName)
    {
        arg_->info().resetName(customName.value());
        return *this;
    }

    ArgCreator<T, TConfig>& operator<<(const ValueName& valueName)
    {
        arg_->info().resetValueName(valueName.value());
        return *this;
    }

    operator T()
    {
        ConfigAccess<TConfig>{cfg_}.addArg(std::move(arg_));
        return T{};
    }

private:
    std::unique_ptr<Arg<T>> arg_;
    TConfig& cfg_;
};

template <typename T, typename TConfig>
ArgCreator<T, TConfig> makeArgCreator(TConfig& cfg,
                                       const std::string& varName,
                                       const std::string& type,
                                       std::function<T&()> argGetter)
{
    return ArgCreator<T, TConfig>{cfg, varName, type, std::move(argGetter)};
}

}

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
#ifndef OPENMC_TALLIES_FILTER_MU_H
#define OPENMC_TALLIES_FILTER_MU_H

#include <vector>

#include <gsl/gsl>

#include "openmc/tallies/filter.h"

namespace openmc {

//==============================================================================
//! Bins the incoming-outgoing direction cosine.  This is only used for scatter
//! reactions.
//==============================================================================

class MuFilter : public Filter
{
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~MuFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override {return "mu";}

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
  const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  void set_bins(gsl::span<double> bins);

private:
  //----------------------------------------------------------------------------
  // Data members

  std::vector<double> bins_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MU_H

max_stars_repo_name: Hit-Weixg/openmc

max_issues_repo_name: mehmeturkmen/openmc

max_forks_repo_name: mehmeturkmen/openmc


text:
/* gsl_histogram2d_oper.c
 * Copyright (C) 2000  Simone Piccardi
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
/***************************************************************
 *
 * File gsl_histogram2d_oper.c:
 * Routine to make operation on 2D histograms.
 * Need GSL library and header.
 * Contains the routines:
 * gsl_histogram2d_same_binning check if two histograms have the same binning
 * gsl_histogram2d_add          add two histogram
 * gsl_histogram2d_sub          subctract two histogram
 * gsl_histogram2d_mult         multiply two histogram
 * gsl_histogram2d_div          divide two histogram
 * gsl_histogram2d_scale        scale histogram contents
 *
 * Author: S. Piccardi
 * Jan. 2000
 *
 ***************************************************************/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram2d.h>

/*
 * gsl_histogram2d_same_binning:
 * control if two histogram have the
 * same binning
 */
int
gsl_histogram2d_equal_bins_p (const gsl_histogram2d * h1,
                              const gsl_histogram2d * h2)
{

  if ((h1->nx != h2->nx) || (h1->ny != h2->ny))
    {
      return 0;
    }
  {
    size_t i;
    /* init ranges */
    for (i = 0; i <= (h1->nx); i++)
      {
        if (h1->xrange[i] != h2->xrange[i])
          {
            return 0;
          }
      }
    for (i = 0; i <= (h1->ny); i++)
      {
        if (h1->yrange[i] != h2->yrange[i])
          {
            return 0;
          }
      }
  }
  return 1;
}

/*
 * gsl_histogram2d_add:
 * add two histogram
 */

int
gsl_histogram2d_add (gsl_histogram2d * h1, const gsl_histogram2d * h2)
{
  size_t i;

  if (!gsl_histogram2d_equal_bins_p (h1, h2))
    {
      GSL_ERROR ("histograms have different binning", GSL_EINVAL);
    }

  for (i = 0; i < (h1->nx) * (h1->ny); i++)
    {
      h1->bin[i] += h2->bin[i];
    }

  return GSL_SUCCESS;
}

/*
 * gsl_histogram2d_sub:
 * subtract two histogram
 */

int
gsl_histogram2d_sub (gsl_histogram2d * h1, const gsl_histogram2d * h2)
{
  size_t i;

  if (!gsl_histogram2d_equal_bins_p (h1, h2))
    {
      GSL_ERROR ("histograms have different binning", GSL_EINVAL);
    }

  for (i = 0; i < (h1->nx) * (h1->ny); i++)
    {
      h1->bin[i] -= h2->bin[i];
    }

  return GSL_SUCCESS;
}

/*
 * gsl_histogram2d_mult:
 * multiply two histogram
 */

int
gsl_histogram2d_mul (gsl_histogram2d * h1, const gsl_histogram2d * h2)
{
  size_t i;

  if (!gsl_histogram2d_equal_bins_p (h1, h2))
    {
      GSL_ERROR ("histograms have different binning", GSL_EINVAL);
    }

  for (i = 0; i < (h1->nx) * (h1->ny); i++)
    {
      h1->bin[i] *= h2->bin[i];
    }

  return GSL_SUCCESS;
}

/*
 * gsl_histogram2d_div:
 * divide two histogram
 */

int
gsl_histogram2d_div (gsl_histogram2d * h1, const gsl_histogram2d * h2)
{
  size_t i;

  if (!gsl_histogram2d_equal_bins_p (h1, h2))
    {
      GSL_ERROR ("histograms have different binning", GSL_EINVAL);
    }

  for (i = 0; i < (h1->nx) * (h1->ny); i++)
    {
      h1->bin[i] /= h2->bin[i];
    }

  return GSL_SUCCESS;
}

/*
 * gsl_histogram2d_scale:
 * scale a histogram by a numeric factor
 */

int
gsl_histogram2d_scale (gsl_histogram2d * h, double scale)
{
  size_t i;

  for (i = 0; i < (h->nx) * (h->ny); i++)
    {
      h->bin[i] *= scale;
    }

  return GSL_SUCCESS;
}

/*
 * gsl_histogram2d_shift:
 * shift a histogram by a numeric offset
 */

int
gsl_histogram2d_shift (gsl_histogram2d * h, double shift)
{
  size_t i;

  for (i = 0; i < (h->nx) * (h->ny); i++)
    {
      h->bin[i] += shift;
    }

  return GSL_SUCCESS;
}


text:
/*
   (c) 2016 Greg Lee
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_sf_bessel.h>

// union cast macro from Matt Gallagher
#define UNION_CAST(x, destType) \
 (((union {__typeof__(x) a; destType b;})x).b)

long int gsl_sf_bessel_J0_c( long int i)
{
   double result = gsl_sf_bessel_J0( UNION_CAST( i, double ) );
   return UNION_CAST( result, long int );
}

/* End of file */

max_stars_repo_name: gleetbone/picolisp-fp

max_issues_repo_name: gleetbone/picolisp-fp

max_forks_repo_name: gleetbone/picolisp-fp


text:
#ifndef SOURCE_UTILS_H_
#define SOURCE_UTILS_H_

// *** source/utils.h ***
// Author: Kevin Wolz, date: 08/2018
//
// Various utilities used in the code.

#include <fstream>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <vector>
#include <string>
#include <unistd.h>
#include <algorithm>
#include <stdlib.h>
#include <sys/stat.h>
#include <time.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_cdf.h>

std::string Timestamp();

double get_signum(double x);

void PrintGslVector(gsl_vector* vector);

void PrintGslVector(gsl_vector* vector, std::string label);

void PrintGslMatrix(gsl_matrix* matrix);

void PrintGslMatrix(gsl_matrix* matrix, std::string label);

bool file_exists(const std::string& filename);

void AppendTwoGslVectors(gsl_vector* vector1, gsl_vector* vector2,
                         gsl_vector* new_vector);

void ConcatenateTwoGslMatrices(gsl_matrix* mmatrix1, gsl_matrix* matrix2,
                               gsl_matrix* new_matrix, bool vertical_axis);

void BlockDiagFromTwoGslMatrices(gsl_matrix* matrix1, gsl_matrix* matrix2,
                                 gsl_matrix* new_matrix, unsigned int dimnew1,
                                 unsigned int dimnew2);

void SubtractTwoDoubleVectors(const std::vector<double> &vector1,
                              const std::vector<double> &vector2,
                              std::vector<double> &difference);

double GslVec1MatVec2(gsl_vector* vector1, gsl_matrix* matrix,
                      gsl_vector* vector2, int dim);

void InvertGslMatrix(gsl_matrix* tobeinverted, unsigned int dim,
                     gsl_matrix* inverted);

double DeterminantGslMatrix(gsl_matrix* square_matrix, unsigned int dim);

double Vec1MatVec2(std::vector<double> vector1, gsl_matrix* matrix,
                   std::vector<double> vector2, int dim);

double TraceGslMatrix(gsl_matrix* square_matrix, int dim);

void MultiplyFourGslMatrices(gsl_matrix* matrix1, gsl_matrix* matrix2,
                             gsl_matrix* matrix3, gsl_matrix* matrix4,
                             int dim, gsl_matrix* resulting_matrix);

void MultiplyThreeGslMatrices(gsl_matrix* matrix1, gsl_matrix* matrix2,
                              gsl_matrix* matrix3, int dim,
                              gsl_matrix* resulting_matrix);

#endif // SOURCE_UTILS_H_

max_stars_repo_name: kevguitar/Robustness

max_issues_repo_name: kevguitar/Robustness

max_forks_repo_name: kevguitar/Robustness


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mersenne.h"
//#include <lapacke/lapacke.h>
#include <stdbool.h>
#include <lapacke.h>
#include <time.h>
#include <sys/time.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif

/* Lattice size and dimensions */
#define ND 4
/* Size of the lattice dimensions */
static int Ldim[ND] = { 4, 4, 4, 4 };

#define NDIRS (2*ND)
int VOLUME;

// Choose update method
#define WORM_UPDATE
//#define LOCAL_UPDATE

//Choose spatial boundary conditions (temporal alway antiperiodic)
//#define PERIODIC
#define ANTIPERIODIC

//Uncomment to include site mass in the determinant
//(does not work with fluctuation determinant)
//#define MASS_IN_MATRIX

/* Flavors */
#define N_FLAVOR 2

/* Choose here the method for calculating the derivative */
//#define FULL_DETERMINANT
#define FLUCTUATION_DETERMINANT

/* Toggle direct measurement for site and link vev */
//#define MEASUREVEV

/* Save memory to fit large lattices
 * Reduces memory use significantly and
 * slightly increases compute times
 */
#define SAVE_MEMORY

/* field values */
#define UNOCCUPIED 0
#define OCCUPIED 1
#define NONE -100   //A meta value for sites that don't exist

/* In Staggered.c */
int site_vector_to_index( int vector[] );
void site_index_to_vector(int index, int vector[] );
void print_config();

/* In fermion_matrix.c */
void Dv( double * y, double * x );
void calc_Dinv( );
double determinant();
void * wrapmalloc( int size, size_t bitsize );

double det_add_monomers(int x1, int x2, int do_flavor[N_FLAVOR]);
double det_move_monomers(int x1, int x2, int do_flavor[N_FLAVOR]);
double det_remove_monomers(int x1, int x2, int do_flavor[N_FLAVOR]);
#ifdef MASS_IN_MATRIX
double det_add_monomers1(int x1, int do_flavor[N_FLAVOR]);
double det_remove_monomers1(int x1, int do_flavor[N_FLAVOR]);
#endif
void update_current_determinant( int do_flavor[N_FLAVOR] );
void update_background( int monomer );

void measure_vev( double * linkvev, double * sitevev );

/* update functions */
double worm_update( int *additions, int *removals, int *m_additions, int *m_removals, int *switches );
void local_update( int *additions, int *removals, int *moves,int *m_additions, int *m_removals, int *m_moves, int *switches );

/* Inversions */
void init_Dinv_cg();
void init_Dinv_deprecated();
void cg( double *source );

max_stars_repo_name: rantahar/Staggered-fourfermion

max_issues_repo_name: rantahar/Staggered-fourfermion

max_forks_repo_name: rantahar/Staggered-fourfermion