text:
//===--- Sudoku/Solver.h                                                ---===//
//
//	Solver functions acting on a Board<Options>&
// Declarations are in Solver.fwd.h
//===----------------------------------------------------------------------===//
#pragma once

#include "Solvers_find.h"
#include "Solvers_remove_option.h"
#include "Solvers_set_option.h"

#include "Sudoku/Board.h"
#include "Sudoku/Location.h"
#include "Sudoku/Location_Utilities.h"
#include "Sudoku/Options.h"
#include "Sudoku/Size.h"
#include "Sudoku/Value.h"
#include <gsl/gsl>

#include <vector>

#include <algorithm>   // minmax
#include <type_traits> // is_base_of

#include "Solver.fwd.h" // Forward declarations

#include <cassert>

namespace Sudoku
{
//	Check if only one option remaining
//	IF true: process answer
template<int N, typename Options>
inline int single_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc)
{
	assert(is_valid(loc));

	if (const Value answer{get_answer(board.at(loc))})
	{
		return single_option(board, loc, answer);
	}
	return 0;
}

//	Process answer from [loc], update board_ options
//	Remove option from rest of row, col and block
template<int N, typename Options>
inline int single_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc,
	const Value value)
{
	{
		assert(is_valid(loc));
		assert(is_valid<N>(value));
		assert(board.at(loc).test(value));
	}
	int changes{};
	changes += set_Value(board, loc, value);
	changes += remove_option_section(board, board.row(loc), loc, value);
	changes += remove_option_section(board, board.col(loc), loc, value);
	changes += remove_option_section(board, board.block(loc), loc, value);
	return changes;
}

// if 2 options in element:
// find exact pair in section:
// remove form other elements in section
template<int N, typename Options>
inline int dual_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc)
{
	using Location = Location<N>;
	{
		assert(is_valid(loc));
		assert(board.at(loc).count() == 2);
	}
	const auto sorted_loc = [loc](const Location L)
	{
		const auto result = std::minmax(loc, L);
		return std::vector<Location>{result.first, result.second};
	};

	int changes{};
	const Options& item{board.at(loc)};
	const auto mask = [](Options x) noexcept
	{
		x[Value{0}] = false;
		return x;
	}(item);

	for (int i{}; i < full_size<N>; ++i)
	{
		// find exact same in board
		if (board[Location(i)] == item && Location(i) != loc)
		{
			// Remove values for rest of shared elements
			if (is_same_row(loc, Location(i)))
			{
				changes += remove_option_section(
					board, board.row(loc), sorted_loc(Location(i)), mask);
			}
			else if (is_same_col(loc, Location(i)))
			{
				changes += remove_option_section(
					board, board.col(loc), sorted_loc(Location(i)), mask);
			}
			// run not if already answered
			if (is_same_block(loc, Location(i)) && !item.is_answer())
			{ // NOTE this is slow
				changes += remove_option_section(
					board, board.block(loc), sorted_loc(Location(i)), mask);
			}
		}
	}
	return changes;
}

// Find subsets in any related sections,
// and remove the values from rest of these sections.
template<int N, typename Options>
inline int multi_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc)
{
	assert(is_valid(loc));

	return multi_option(board, loc, board.at(loc).count());
}

// Find subsets in any related section,
// and remove the values from rest of these sections.
template<int N, typename Options>
constexpr int multi_option(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc,
	const size_t count)
{
	assert(is_valid(loc));

	// set execution domain & specializations
	switch (count)
	{
	case 0: [[fallthrough]];     // already answered
	case elem_size<N>: return 0; // no result anyway
	case 1: return single_option(board, loc);
	case 2: return dual_option(board, loc);
	default: assert(count < elem_size<N>); break;
	}

	int changes{};                      // performance counter
	const Options& item{board.at(loc)}; // input item, to match with
	assert(item.count() == count);
	const auto mask = [](Options x) noexcept
	{
		x[Value{0}] = false;
		return x;
	}(item);

	const auto list = list_where_subset(board, item);

	// Further select per section.
	// When (subset size == #options) then: all answers in this selection.
	// Therefore: Remove values for rest of section.
	if (const auto in_row{get_same_row(loc, list)}; in_row.size() == count)
	{
		changes += remove_option_section(board, board.row(loc), in_row, mask);
	}
	if (const auto in_col{get_same_col(loc, list)}; in_col.size() == count)
	{
		changes += remove_option_section(board, board.col(loc), in_col, mask);
	}
	if (const auto in_block{get_same_block(loc, list)};
		in_block.size() == count)
	{
		// NOTE this is slow
		changes +=
			remove_option_section(board, board.block(loc), in_block, mask);
	}
	return changes;
}

// Solver: Find options appearing only once in a section and set as answer
template<int N, typename Options, typename SectionT>
inline int unique_in_section(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
	}
	const auto worker = appearance_once<N>(section);
	return set_uniques(board, section, worker);
}

//	Solver: find and process values appearing 1 to base_size times in section:
//	[row/col] IF all in same block -> remove from rest of block
//	[block] IF all in same row/col -> remove from rest of row/col
template<int N, typename Options, typename SectionT>
inline int section_exclusive(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<typename SectionT::value_type, Options>);
	}
	int changes{}; // performance counter

	size_t i{2};
	auto appearing             = appearance_sets<N>(section);
	const auto renew_appearing = [&i, &a = appearing, &section]()
	{
		i = 2;
		a = appearance_sets<N>(section);
	};
	while (i < appearing.size()) // won't run if condition fails
	{
		// unique specialization
		if (appearing[1].count_all() > 0)
		{
			changes += set_uniques(board, section, appearing[1]);
			renew_appearing();
			continue;
		}
		if (appearing.at(i).count_all() > 0)
		{ // for [row/col]: if in same block: remove from rest block
			// for [block]: if in same row/col: remove from rest row/col
			if (const int tmp_ = set_section_locals(
					board, section, gsl::narrow_cast<int>(i), appearing.at(i)))
			{
				changes += tmp_;
				renew_appearing();
				continue;
			}
		}
		++i;
	}
	return changes;
}

// cell containing 2 options
//	find in board_	size() == 2
//		find exact same in rest board
//			IF pair in same row/col -> remove values from rest row/col
//			IF same block -> remove values from rest block
//		repeat until end of board_
//	repeat until end of board_

// cell containing 2 options
//	find in board_	size() == 2
//		find exact same in rest row
//			IF found:	remove values from rest row
//		find exact same in rest col
//		find exact same in rest block

} // namespace Sudoku

max_stars_repo_name: Farwaykorse/fwkSudoku

max_issues_repo_name: Farwaykorse/fwkSudoku

max_forks_repo_name: Farwaykorse/fwkSudoku


text:
/* linalg/lu.c
 *
 * Copyright (C) 2020 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.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

static int LU_band_decomp_L2 (const size_t M, const size_t lb, const size_t ub,
                              gsl_matrix * AB, gsl_vector_uint * ipiv);

/* Factorise a (p,q) banded M x N matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.
 *
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * Inputs: M   - number of rows in matrix
 *         lb  - lower bandwidth
 *         ub  - upper bandwidth
 *         AB  - matrix in band storage format, N-by-(2*p + q + 1)
 *         piv - pivot vector, size MIN(M, N)
 */

int
gsl_linalg_LU_band_decomp (const size_t M, const size_t lb, const size_t ub, gsl_matrix * AB, gsl_vector_uint * piv)
{
  const size_t N = AB->size1;
  const size_t minMN = GSL_MIN(M, N);

  if (lb >= M)
    {
      GSL_ERROR ("lower bandwidth must be less than M", GSL_EDOM);
    }
  else if (ub >= N)
    {
      GSL_ERROR ("upper bandwidth must be less than N", GSL_EDOM);
    }
  else if (AB->size2 != 2*lb + ub + 1)
    {
      GSL_ERROR ("matrix size inconsistent with bandwidths", GSL_EBADLEN);
    }
  else if (piv->size != minMN)
    {
      GSL_ERROR ("pivot vector must have length MIN(M,N)", GSL_EBADLEN);
    }
  else
    {
      int status;

      status = LU_band_decomp_L2 (M, lb, ub, AB, piv);

      return status;
    }
}

int
gsl_linalg_LU_band_solve (const size_t lb, const size_t ub, const gsl_matrix * LUB,
                          const gsl_vector_uint * piv, const gsl_vector * b, gsl_vector * x)
{
  const size_t N = LUB->size1;

  if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (N != b->size)
    {
      GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
    }
  else if (lb >= N)
    {
      GSL_ERROR ("lower bandwidth must be less than N", GSL_EDOM);
    }
  else if (ub >= N)
    {
      GSL_ERROR ("upper bandwidth must be less than N", GSL_EDOM);
    }
  else if (LUB->size2 != 2*lb + ub + 1)
    {
      GSL_ERROR ("matrix size inconsistent with bandwidths", GSL_EBADLEN);
    }
  else if (piv->size != N)
    {
      GSL_ERROR ("pivot vector must have length N", GSL_EBADLEN);
    }
  else
    {
      int status;

      gsl_vector_memcpy(x, b);

      status = gsl_linalg_LU_band_svx(lb, ub, LUB, piv, x);

      return status;
    }
}

int
gsl_linalg_LU_band_svx (const size_t lb, const size_t ub, const gsl_matrix * LUB,
                        const gsl_vector_uint * piv, gsl_vector * x)
{
  const size_t N = LUB->size1;

  if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else if (lb >= N)
    {
      GSL_ERROR ("lower bandwidth must be less than N", GSL_EDOM);
    }
  else if (ub >= N)
    {
      GSL_ERROR ("upper bandwidth must be less than N", GSL_EDOM);
    }
  else if (LUB->size2 != 2*lb + ub + 1)
    {
      GSL_ERROR ("matrix size inconsistent with bandwidths", GSL_EBADLEN);
    }
  else if (piv->size != N)
    {
      GSL_ERROR ("pivot vector must have length N", GSL_EBADLEN);
    }
  else
    {
      if (lb > 0)
        {
          size_t j;

          for (j = 0; j < N - 1; ++j)
            {
              size_t pj = gsl_vector_uint_get(piv, j);
              double * xj = gsl_vector_ptr(x, j);
              size_t lm = GSL_MIN(lb, N - j - 1);
              gsl_vector_view xv = gsl_vector_subvector(x, j + 1, lm);
              gsl_vector_const_view yv = gsl_matrix_const_subrow(LUB, j, lb + ub + 1, lm);

              if (j != pj)
                {
                  double xl = gsl_vector_get(x, pj);
                  gsl_vector_set(x, pj, *xj);
                  *xj = xl;
                }

              gsl_blas_daxpy(-(*xj), &yv.vector, &xv.vector);
            }
        }

      /* solve U x = b */
      cblas_dtbsv(CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit,
                  (int) N, (int) (lb + ub), LUB->data, LUB->tda,
                  x->data, x->stride);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_band_unpack (const size_t M, const size_t lb, const size_t ub, const gsl_matrix * LUB,
                           const gsl_vector_uint * piv, gsl_matrix * L, gsl_matrix * U)
{
  const size_t N = LUB->size1;
  const size_t minMN = GSL_MIN(M, N);

  if (ub >= N)
    {
      GSL_ERROR ("upper bandwidth must be < N", GSL_EDOM);
    }
  else if (lb >= M)
    {
      GSL_ERROR ("lower bandwidth must be < M", GSL_EDOM);
    }
  else if (LUB->size2 != 2*lb + ub + 1)
    {
      GSL_ERROR ("matrix size inconsistent with bandwidths", GSL_EBADLEN);
    }
  else if (piv->size != minMN)
    {
      GSL_ERROR ("pivot vector must have length MIN(M,N)", GSL_EBADLEN);
    }
  else if (L->size1 != M || L->size2 != minMN)
    {
      GSL_ERROR ("L matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (U->size1 != minMN || U->size2 != N)
    {
      GSL_ERROR ("U matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      const size_t ub_U = lb + ub;
      size_t j;

      gsl_matrix_set_identity(L);
      gsl_matrix_set_zero(U);

      /* compute L */
      if (lb > 0)
        {
          const size_t jstart = (M > N) ? minMN : minMN - 1;
          size_t j;

          for (j = jstart; j > 0 && j--; )
            {
              size_t pj = gsl_vector_uint_get(piv, j);
              size_t lm = GSL_MIN(lb, M - j - 1);
              gsl_vector_const_view xv = gsl_matrix_const_subrow(LUB, j, lb + ub + 1, lm);
              gsl_vector_const_view yv = gsl_matrix_const_subrow(L, j, 0, minMN);
              gsl_matrix_view m = gsl_matrix_submatrix(L, j + 1, 0, lm, minMN);

              gsl_blas_dger(1.0, &xv.vector, &yv.vector, &m.matrix);

              if (j != pj)
                {
                  gsl_vector_view Lj = gsl_matrix_row(L, j);
                  gsl_vector_view Lpj = gsl_matrix_row(L, pj);
                  gsl_blas_dswap(&Lj.vector, &Lpj.vector);
                }
            }
        }

      /* fill in U */
      for (j = 0; j <= GSL_MIN(ub_U, N - 1); ++j)
        {
          gsl_vector_const_view src = gsl_matrix_const_subcolumn(LUB, ub_U - j, j, GSL_MIN(M, N - j));
          gsl_vector_view dest = gsl_matrix_superdiagonal(U, j);
          gsl_vector_memcpy(&dest.vector, &src.vector);
        }

      return GSL_SUCCESS;
    }
}

/*
LU_band_decomp_L2
  LU decomposition with partial pivoting using Level 2 BLAS

Inputs: M    - number of rows in matrix
        lb   - lower bandwidth
        ub   - upper bandwidth
        AB   - on input, matrix to be factored; on output, L and U factors
               N-by-(2*p + q + 1)
        ipiv - (output) array containing row swaps

Notes:
1) Based on LAPACK DGBTF2
*/

static int
LU_band_decomp_L2 (const size_t M, const size_t lb, const size_t ub,
                   gsl_matrix * AB, gsl_vector_uint * ipiv)
{
  const size_t N = AB->size1;
  const size_t minMN = GSL_MIN(M, N);

  if (ipiv->size != minMN)
    {
      GSL_ERROR ("ipiv length must equal MIN(M,N)", GSL_EBADLEN);
    }
  else if (lb >= M)
    {
      GSL_ERROR ("lower bandwidth must be less than M", GSL_EDOM);
    }
  else if (ub >= N)
    {
      GSL_ERROR ("upper bandwidth must be less than N", GSL_EDOM);
    }
  else if (AB->size2 != 2*lb + ub + 1)
    {
      GSL_ERROR ("matrix size inconsistent with bandwidths", GSL_EBADLEN);
    }
  else
    {
      const size_t ub_U = lb + ub; /* upper bandwidth of U factor */
      const size_t ldab = AB->size2;
      size_t ju = 0;
      size_t j;

      if (lb > 0)
        {
          /* initialize fill-in elements to zero */
          gsl_matrix_view m = gsl_matrix_submatrix(AB, 0, 0, N, lb);
          gsl_matrix_set_zero(&m.matrix);
        }

      for (j = 0; j < minMN; ++j)
        {
          size_t lbj = GSL_MIN(lb, M - j - 1); /* subdiagonal elements in column j */
          gsl_vector_view x = gsl_matrix_subrow(AB, j, ub_U, lbj + 1);
          gsl_vector_view y;
          CBLAS_INDEX_t j_pivot = gsl_blas_idamax(&x.vector);
          double * ptr;

          gsl_vector_uint_set(ipiv, j, j + j_pivot);

          ptr = gsl_matrix_ptr(AB, j, ub_U + j_pivot);
          if (*ptr != 0.0)
            ju = GSL_MAX(ju, GSL_MIN(j + ub + j_pivot, N - 1));

          if (j_pivot != 0)
            {
              double *ptr2;

              /* swap columns */

              x = gsl_vector_view_array_with_stride(ptr, ldab - 1, ju - j + 1);

              ptr2 = gsl_matrix_ptr(AB, j, ub_U);
              y = gsl_vector_view_array_with_stride(ptr2, ldab - 1, ju - j + 1);

              gsl_blas_dswap(&x.vector, &y.vector);
            }

          if (lbj > 0)
            {
              double tmp = gsl_matrix_get(AB, j, ub_U);

              x = gsl_matrix_subrow(AB, j, ub_U + 1, lbj);
              gsl_blas_dscal(1.0 / tmp, &x.vector);

              if (ju > j)
                {
                  gsl_matrix_view m = gsl_matrix_submatrix(AB, j + 1, ub_U, ju - j, lbj);

                  ptr = gsl_matrix_ptr(AB, j + 1, ub_U - 1);
                  y = gsl_vector_view_array_with_stride(ptr, ldab - 1, ju - j);

                  m.matrix.tda = ldab - 1; /* unskew matrix */
                  gsl_blas_dger(-1.0, &y.vector, &x.vector, &m.matrix);
                }
            }
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* linalg/ldlt.c
 *
 * Copyright (C) 2018 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.
 */

/* L D L^T decomposition of a symmetric positive semi-definite matrix */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

static double ldlt_norm1(const gsl_matrix * A);
static int ldlt_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);

/*
gsl_linalg_ldlt_decomp()
  Perform L D L^T decomposition of a symmetric positive
semi-definite matrix using lower triangle

Inputs: A - (input) symmetric, positive semi-definite matrix
            (output) lower triangle contains L factor;
                     diagonal contains D

Return: success/error

Notes:
1) Based on algorithm 4.1.1 of Golub and Van Loan, Matrix Computations (4th ed).
2) The first subrow A(1, 2:end) is used as temporary workspace
3) The 1-norm ||A||_1 of the original matrix is stored in the upper right corner on output
*/

int
gsl_linalg_ldlt_decomp (gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("LDLT decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t i, j;
      double a00, anorm;
      gsl_vector_view work, v;

      /* check for quick return */
      if (N == 1)
        return GSL_SUCCESS;

      /* compute ||A||_1 */
      anorm = ldlt_norm1(A);

      /* special case first column */
      a00 = gsl_matrix_get(A, 0, 0);
      if (a00 == 0.0)
        {
          GSL_ERROR ("matrix is singular", GSL_EDOM);
        }

      v = gsl_matrix_subcolumn(A, 0, 1, N - 1);
      gsl_vector_scale(&v.vector, 1.0 / a00);

      /* use first subrow A(1, 2:end) as temporary workspace */
      work = gsl_matrix_subrow(A, 0, 1, N - 1);

      for (j = 1; j < N; ++j)
        {
          gsl_vector_view w = gsl_vector_subvector(&work.vector, 0, j);
          double ajj = gsl_matrix_get(A, j, j);
          double dval;

          for (i = 0; i < j; ++i)
            {
              double aii = gsl_matrix_get(A, i, i);
              double aji = gsl_matrix_get(A, j, i);
              gsl_vector_set(&w.vector, i, aji * aii);
            }

          v = gsl_matrix_subrow(A, j, 0, j); /* A(j,1:j-1) */
          gsl_blas_ddot(&v.vector, &w.vector, &dval);
          ajj -= dval;

          if (ajj == 0.0)
            {
              GSL_ERROR ("matrix is singular", GSL_EDOM);
            }

          gsl_matrix_set(A, j, j, ajj);

          if (j < N - 1)
            {
              double ajjinv = 1.0 / ajj;
              gsl_matrix_view m = gsl_matrix_submatrix(A, j + 1, 0, N - j - 1, j); /* A(j+1:n, 1:j-1) */
              v = gsl_matrix_subcolumn(A, j, j + 1, N - j - 1);                    /* A(j+1:n, j) */
              gsl_blas_dgemv(CblasNoTrans, -ajjinv, &m.matrix, &w.vector, ajjinv, &v.vector);
            }
        }

      /* save ||A||_1 in upper right corner */
      gsl_matrix_set(A, 0, N - 1, anorm);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_ldlt_solve (const gsl_matrix * LDLT,
                       const gsl_vector * b,
                       gsl_vector * x)
{
  if (LDLT->size1 != LDLT->size2)
    {
      GSL_ERROR ("LDLT matrix must be square", GSL_ENOTSQR);
    }
  else if (LDLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LDLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* copy x <- b */
      gsl_vector_memcpy (x, b);

      status = gsl_linalg_ldlt_svx(LDLT, x);

      return status;
    }
}

int
gsl_linalg_ldlt_svx (const gsl_matrix * LDLT,
                     gsl_vector * x)
{
  if (LDLT->size1 != LDLT->size2)
    {
      GSL_ERROR ("LDLT matrix must be square", GSL_ENOTSQR);
    }
  else if (LDLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_const_view diag = gsl_matrix_const_diagonal(LDLT);

      /* solve for z using forward-substitution, L z = b */
      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasUnit, LDLT, x);

      /* solve for y, D y = z */
      gsl_vector_div(x, &diag.vector);

      /* perform back-substitution, L^T x = y */
      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasUnit, LDLT, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_ldlt_rcond (const gsl_matrix * LDLT, double * rcond,
                       gsl_vector * work)
{
  const size_t M = LDLT->size1;
  const size_t N = LDLT->size2;

  if (M != N)
    {
      GSL_ERROR ("LDLT matrix must be square", GSL_ENOTSQR);
    }
  else if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      int status;
      double Anorm;    /* ||A||_1 */
      double Ainvnorm; /* ||A^{-1}||_1 */

      if (N == 1)
        Anorm = fabs(gsl_matrix_get(LDLT, 0, 0));
      else
        Anorm = gsl_matrix_get(LDLT, 0, N - 1);

      *rcond = 0.0;

      /* don't continue if matrix is singular */
      if (Anorm == 0.0)
        return GSL_SUCCESS;

      /* estimate ||A^{-1}||_1 */
      status = gsl_linalg_invnorm1(N, ldlt_Ainv, (void *) LDLT, &Ainvnorm, work);

      if (status)
        return status;

      if (Ainvnorm != 0.0)
        *rcond = (1.0 / Anorm) / Ainvnorm;

      return GSL_SUCCESS;
    }
}

/* compute 1-norm of symmetric matrix stored in lower triangle */
static double
ldlt_norm1(const gsl_matrix * A)
{
  const size_t N = A->size1;
  double max = 0.0;
  size_t i, j;

  for (j = 0; j < N; ++j)
    {
      gsl_vector_const_view v = gsl_matrix_const_subcolumn(A, j, j, N - j);
      double sum = gsl_blas_dasum(&v.vector);

      /* now add symmetric elements from above diagonal */
      for (i = 0; i < j; ++i)
        {
          double Aij = gsl_matrix_get(A, i, j);
          sum += fabs(Aij);
        }

      if (sum > max)
        max = sum;
    }

  return max;
}

/* x := A^{-1} x = A^{-t} x, A = L D L^T */
static int
ldlt_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
{
  int status;
  gsl_matrix * A = (gsl_matrix * ) params;
  gsl_vector_const_view diag = gsl_matrix_const_diagonal(A);

  (void) TransA; /* unused parameter warning */

  /* compute L^{-1} x */
  status = gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasUnit, A, x);
  if (status)
    return status;

  /* compute D^{-1} x */
  gsl_vector_div(x, &diag.vector);

  /* compute L^{-t} x */
  status = gsl_blas_dtrsv(CblasLower, CblasTrans, CblasUnit, A, x);
  if (status)
    return status;

  return GSL_SUCCESS;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * 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 SIM_PHYSICS_DVL_AERO_DATABASE_H_
#define SIM_PHYSICS_DVL_AERO_DATABASE_H_

#include <gsl/gsl_vector.h>
#include <stddef.h>
#include <stdint.h>

#include <string>
#include <vector>

#include "common/c_math/force_moment.h"
#include "common/c_math/linalg.h"
#include "common/c_math/vec3.h"
#include "common/macros.h"
#include "sim/physics/aero_types.h"
#include "system/labels.h"

class DvlAeroDatabase {
 public:
  explicit DvlAeroDatabase(const std::string &filename);
  ~DvlAeroDatabase();

  // Returns the total force-moment coefficient at a specific (alpha,
  // beta, omega_hat, flap deflection) operating point.
  void CalcForceMomentCoeff(double alpha, double beta, const Vec3 &omega_hat,
                            const Vec &flaps, ForceMoment *cfm,
                            double thrust_coeff) const {
    CalcForceMomentCoeff(alpha, beta, omega_hat, flaps, cfm, thrust_coeff,
                         nullptr);
  }

  // See function above.  Optionally returns the breakdown of the
  // force-moment coefficient into control and stability derivatives
  // if aero_coeffs is not null.
  void CalcForceMomentCoeff(double alpha, double beta, const Vec3 &omega_hat,
                            const Vec &flaps, ForceMoment *cfm,
                            double thrust_coeff,
                            DvlAeroCoeffs *aero_coeffs) const;

  int32_t num_flaps() const { return num_flaps_; }
  double reynolds_number() const { return reynolds_number_; }

 private:
  // Returns true if the database passes basic sign checks.
  bool IsValid() const;

  // Returns the total force-moment coefficient at a specific (alpha,
  // beta, omega_hat, flap deflection) operating point and the
  // breakdown of this value into control and stability derivatives.
  void CalcDvlAeroCoeffs(double alpha, double beta, const Vec &flaps,
                         DvlAeroCoeffs *coeffs) const;

  void LookupForceMomentCoeff2D(const gsl_vector &cfm_data, size_t size1,
                                size_t size2, int32_t s0, int32_t t0, double ds,
                                double dt, ForceMoment *cfm) const;
  void LookupForceMomentCoeff3D(const gsl_vector &cfm_data, int32_t indices[],
                                double du, double dv, double dw,
                                ForceMoment *cfm) const;

  // Number of flaps used in database.
  int32_t num_flaps_;

  // Reynolds number [#] at which the database was calculated.
  double reynolds_number_;

  // List of angles-of-attack [rad] used in database look-up table.
  gsl_vector *alphas_;

  // List of sideslip angles [rad] used in database look-up table.
  gsl_vector *betas_;

  // Vector of lists of flap deflections [rad], one list for each
  // flap, used in database look-up table.
  std::vector<gsl_vector *> deltas_;

  // Tensor grid of force-moment coefficients at zero flap deflections
  // over angle-of-attack and sideslip, stored as a vector.
  gsl_vector *cfm_;

  // Tensor grid of force-moment stability derivatives at zero flap
  // deflections over angle-of-attack and sideslip, stored as a
  // vector.
  gsl_vector *dcfm_dp_;
  gsl_vector *dcfm_dq_;
  gsl_vector *dcfm_dr_;

  // Vector of tensor grids of increments to the force-moment
  // coefficients, due to flap deflections, for each flap.  The tensor
  // grid is over angle-of-attack, sideslip, and flap deflection, and
  // is stored as a vector.
  std::vector<gsl_vector *> flap_cfms_;

  // Vector of tensor grids of increments to the force-moment
  // stability derivatives, due to flap deflections, for each flap.
  // The tensor grid is over angle-of-attack, sideslip, and flap
  // deflection, and is stored as a vector.
  std::vector<gsl_vector *> flap_dcfm_dps_;
  std::vector<gsl_vector *> flap_dcfm_dqs_;
  std::vector<gsl_vector *> flap_dcfm_drs_;

  DISALLOW_COPY_AND_ASSIGN(DvlAeroDatabase);
};

#endif  // SIM_PHYSICS_DVL_AERO_DATABASE_H_

ext: h


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef rtree_77e9ff20_6040_46f3_b87d_ab64a4ac506a_h
#define rtree_77e9ff20_6040_46f3_b87d_ab64a4ac506a_h

#include <gslib/error.h>
#include <gslib/tree.h>
#include <gslib/utility.h>

__gslib_begin__

class __gs_novtable rtree_spec abstract
{
public:
    enum type
    {
        spec_line,
        /* more... */
    };

    virtual ~rtree_spec() {}
    virtual type get_type() const = 0;
    virtual bool query_overlap(const rectf& rc) const = 0;
    virtual bool query_overlap(const pointf& p1, const pointf& p2) const = 0;
};

class rtree_spec_line:
    public rtree_spec
{
public:
    rtree_spec_line(const pointf& p1, const pointf& p2): _p1(p1), _p2(p2) {}
    virtual type get_type() const override { return spec_line; }
    virtual bool query_overlap(const rectf& rc) const override { return is_line_rect_overlapped(_p1, _p2, rc); }
    virtual bool query_overlap(const pointf& p1, const pointf& p2) const override
    {
        vec2 intsp, d1, d2;
        d1.sub(_p2, _p1);
        d2.sub(p2, p1);
        if(is_parallel(d1, d2))
            return false;
        intersectp_linear_linear(intsp, _p1, p1, d1, d2);
        float t1 = linear_reparameterize(_p1, _p2, intsp);
        if(t1 < 0.f || t1 > 1.f)
            return false;
        float t2 = linear_reparameterize(p1, p2, intsp);
        return t2 >= 0.f && t2 <= 1.f;
    }

private:
    pointf          _p1, _p2;
};

template<class _bind, int _bsize = sizeof(_bind)>
class rtree_entity
{
public:
    typedef rtree_entity<_bind, _bsize> myref;
    typedef _bind bind_type;
    typedef const _bind& bind_arg;

protected:
    rectf           _rect;
    bind_type       _bind_arg;
    rtree_spec*     _spec = nullptr;

public:
    rtree_entity() {}
    rtree_entity(bind_arg ba, const rectf& rc): _rect(rc), _bind_arg(ba) {}
    ~rtree_entity() { set_spec(nullptr); }
    const rectf& const_rect() const { return _rect; }
    rectf& get_rect() { return _rect; }
    bind_arg get_bind_arg() const { return _bind_arg; }
    rtree_spec* get_spec() const { return _spec; }
    void set_rect(const rectf& rc) { _rect = rc; }
    void set_bind_arg(bind_arg ba) { _bind_arg = ba; }
    void set_spec(rtree_spec* spec)
    {
        if(_spec)
            delete _spec;
        _spec = spec;
    }
    rtree_spec* detach_spec()
    {
        auto* spec = _spec;
        _spec = nullptr;
        return spec;
    }
};

template<class _bind>
class rtree_entity<_bind, sizeof(void*)>
{
public:
    typedef rtree_entity<_bind, sizeof(void*)> myref;
    typedef _bind bind_type;
    typedef _bind bind_arg;

protected:
    rectf           _rect;
    bind_type       _bind_arg = 0;
    rtree_spec*     _spec = nullptr;

public:
    rtree_entity() {}
    rtree_entity(bind_arg ba, const rectf& rc): _rect(rc), _bind_arg(ba) {}
    ~rtree_entity() { set_spec(nullptr); }
    const rectf& const_rect() const { return _rect; }
    rectf& get_rect() { return _rect; }
    bind_arg get_bind_arg() const { return _bind_arg; }
    rtree_spec* get_spec() const { return _spec; }
    void set_rect(const rectf& rc) { _rect = rc; }
    void set_bind_arg(bind_arg ba) { _bind_arg = ba; }
    void set_spec(rtree_spec* spec)
    {
        if(_spec)
            delete _spec;
        _spec = spec;
    }
    rtree_spec* detach_spec()
    {
        auto* spec = _spec;
        _spec = nullptr;
        return spec;
    }
};

template<class _entity>
using rtree_node = _treenode_cpy_wrapper<_entity>;

template<class _rect>
static float calc_enlarge_area(const _rect& rc1, const _rect& rc2, float s)
{
    assert(s == rc1.area());
    float s2 = rc2.area();
    _rect u, d;
    union_rect(u, rc1, rc2);
    if(intersect_rect(d, rc1, rc2))
        return u.area() - s - s2 + d.area();
    return u.area() - s - s2;
}

template<class _rect>
static float calc_enlarge_area(const _rect& rc1, const _rect& rc2)
{
    return calc_enlarge_area(rc1, rc2, rc1.area());
}

inline float calc_orientation(const float coef[3], const pointf& p)
{
    assert(coef);
    return coef[0] * p.x + coef[1] * p.y + coef[2];
}

inline void calc_perpendicular_bisector(float coef[3], const pointf& p1, const pointf& p2)
{
    assert(coef);
    pointf c;
    c.x = (p1.x + p2.x) * 0.5f;
    c.y = (p1.y + p2.y) * 0.5f;
    float d[2];
    d[0] = p1.y - p2.y;
    d[1] = p1.x - p2.x;
    if(abs(d[0]) < 1e-5f) {
        coef[0] = 1.f;
        coef[1] = -c.x;
        coef[2] = 0.f;
    }
    else if(abs(d[1]) < 1e-5f) {
        coef[0] = -c.y;
        coef[1] = 1.f;
        coef[2] = 0.f;
    }
    else {
        coef[0] = d[1] / d[0];
        coef[1] = 1.f;
        coef[2] = -c.x * coef[0] - c.y;
    }
    /* test orientation */
    float t = calc_orientation(coef, p1);
    if(t > 0.f) {
        coef[0] = -coef[0];
        coef[1] = -coef[1];
        coef[2] = -coef[2];
    }
}

template<class _tree>
class rtree_finder
{
public:
    typedef _tree mytree;
    typedef typename _tree::value value;
    typedef typename _tree::iterator iterator;
    typedef typename _tree::const_iterator const_iterator;
    typedef typename value::bind_arg bind_arg;

protected:
    mytree&         _mytree;

public:
    rtree_finder(mytree& mt): _mytree(mt) {}
    iterator find(bind_arg ba, const rectf& rc) const { return find(_mytree.get_root(), ba, rc); }

protected:
    iterator find(iterator p, bind_arg ba, const rectf& rc) const
    {
        assert(p.is_valid());
        if(!is_rect_intersected(rc, p->const_rect()))
            return iterator(nullptr);
        if(p.is_leaf())
            return (p->get_bind_arg() == ba) ? p : iterator(nullptr);
        for(auto i = p.child(); i.is_valid(); i.to_next()) {
            auto f = find(i, ba, rc);
            if(f.is_valid())
                return f;
        }
        return iterator(nullptr);
    }
};

/*
 * Here we implement an algorithm for the bulk loading. For a set of input inserting into
 * the rtree, original insert was a waste.
 * The algorithm was very similar to the STR insert as the following paper stated:
 * http://www.dtic.mil/dtic/tr/fulltext/u2/a324493.pdf
 */
template<int _max_record, int _min_record, class _tree>
class rtree_str_load
{
public:
    typedef _tree mytree;
    typedef typename _tree::value value;
    typedef typename _tree::wrapper wrapper;
    typedef typename _tree::iterator iterator;
    typedef typename _tree::const_iterator const_iterator;
    typedef typename _tree::value value;
    typedef typename value::bind_arg bind_arg;
    typedef list<value> value_list;
    static const int max_record = _max_record;
    static const int min_record = _min_record;
    struct bulk_node;
    typedef list<bulk_node*> bulk_list;
    typedef list<wrapper*> wrapper_list;

    struct bulk_node:
        public value
    {
        pointf      center;
        wrapper*    assoc_wrapper;

    public:
        bulk_node(const value& v):
            value(v)
        {
            auto& rc = v.const_rect();
            center.x = 0.5f * (rc.left + rc.right);
            center.y = 0.5f * (rc.top + rc.bottom);
            assoc_wrapper = nullptr;
        }
        bulk_node(wrapper* w):
            value(w->get_ref())
        {
            auto& v = w->get_ref();
            auto& rc = v.const_rect();
            center.x = 0.5f * (rc.left + rc.right);
            center.y = 0.5f * (rc.top + rc.bottom);
            assoc_wrapper = w;
        }
    };

public:
    rtree_str_load(mytree& tr): _mytree(tr) {}
    template<class _vsl>
    void load(const _vsl& input)
    {
        wrapper_list in, out;
        for(auto& p : input) {
            auto* w = new wrapper;
            w->born();
            auto& v = w->get_ref();
            v.set_bind_arg(p.get_bind_arg());
            v.set_rect(p.const_rect());
            v.set_spec(p.detach_spec());
            in.push_back(w);
        }
        for(;;) {
            int s = packing(out, in);
            if(s == 1) {
                _mytree.set_root(out.front());
                return;
            }
            out.swap(in);
            out.clear();
        }
    }

protected:
    mytree&         _mytree;

protected:
    void set_as_children(wrapper* p, wrapper_list& ch)
    {
        assert(p);
        assert(!ch.empty());
        int size = (int)ch.size();
        if(size == 1) {
            p->acquire_children(wrapper::children(ch.front()));
            return;
        }
        auto i = ch.begin(), j = std::next(i), end = ch.end();
        for(; j != end; i = j ++)
            wrapper::children::join(*i, *j);
        p->acquire_children(wrapper::children(ch.front(), ch.back(), (int)ch.size()));
    }
    void set_bound_rect(wrapper* p)
    {
        assert(p);
        float l = FLT_MAX, t = FLT_MAX, r = -FLT_MAX, b = -FLT_MAX;
        p->for_each([&](wrapper* w) {
            assert(w);
            auto& v = w->get_ref();
            auto& rc = v.get_rect();
            l = gs_min(l, rc.left);
            t = gs_min(t, rc.top);
            r = gs_max(r, rc.right);
            b = gs_max(b, rc.bottom);
        });
        p->born();
        auto& v = p->get_ref();
        rectf rc;
        rc.set_ltrb(l, t, r, b);
        v.set_rect(rc);
    }
    int packing(wrapper_list& output, wrapper_list& input)
    {
        assert(output.empty());
        assert(!input.empty());
        int size = (int)input.size();
        if(size <= max_record) {
            auto* p = new wrapper;
            set_as_children(p, input);
            set_bound_rect(p);
            output.push_back(p);
            return 1;
        }
        list<bulk_node> bulks;
        for(auto* w : input)
            bulks.push_back(bulk_node(w));
        bulks.sort([](const bulk_node& n1, const bulk_node& n2)-> bool {
            return n1.center.x < n2.center.x;
        });
        float rough_slice_count = sqrtf((float)size / max_record);
        float rough_slice_size = (float)size / rough_slice_count;
        int slice_size = (int)(rough_slice_size + 0.5f);
        int full_slice_count = size / slice_size;
        int slice_count = (size % slice_size) ? full_slice_count + 1 : full_slice_count;
        auto* slices = new list<bulk_node*>[slice_count];
        assert(slices);
        auto i = bulks.begin(), end = bulks.end();
        bool ascend = true;
        auto asc_pr = [](const bulk_node* n1, const bulk_node* n2)-> bool { return n1->center.y < n2->center.y; };
        auto desc_pr = [](const bulk_node* n1, const bulk_node* n2)-> bool { return n1->center.y > n2->center.y; };
        for(int j = 0, k = 0; i != end; ++ i) {
            assert(k < slice_count);
            slices[k].push_back(&(*i));
            if(++ j >= slice_size) {
                assert((int)slices[k].size() == slice_size);
                ascend ? slices[k].sort(asc_pr) : slices[k].sort(desc_pr);
                ascend = !ascend;
                j = 0, k ++;
            }
        }
        if(slice_count != full_slice_count) {
            assert(full_slice_count + 1 == slice_count);
            ascend ? slices[full_slice_count].sort(asc_pr) : slices[full_slice_count].sort(desc_pr);
        }
        /* join the list */
        vector<bulk_node*> sorted;
        for(int j = 0; j < slice_count; j ++)
            sorted.insert(sorted.end(), slices[j].begin(), slices[j].end());
        divide(output, sorted);
        delete [] slices;
        return (int)output.size();
    }
    void divide(wrapper_list& output, vector<bulk_node*>& sorted)
    {
        int size = (int)sorted.size();
        for(int i = 0; i < size; i += max_record) {
            int left = size - i;
            if(left < max_record + min_record) {
                if(left <= max_record) {
                    wrapper_list ch;
                    for(int j = i; j < size; j ++)
                        ch.push_back(sorted.at(j)->assoc_wrapper);
                    auto* w = new wrapper;
                    set_as_children(w, ch);
                    set_bound_rect(w);
                    output.push_back(w);
                }
                else {
                    int size1 = left / 2;
                    int size2 = left - size1;
                    wrapper_list ch1, ch2;
                    for(int j = 0; j < size1; j ++)
                        ch1.push_back(sorted.at(j + i)->assoc_wrapper);
                    for(int j = i + size1; j < size; j ++)
                        ch2.push_back(sorted.at(j)->assoc_wrapper);
                    auto* w1 = new wrapper;
                    auto* w2 = new wrapper;
                    set_as_children(w1, ch1);
                    set_bound_rect(w1);
                    set_as_children(w2, ch2);
                    set_bound_rect(w2);
                    output.push_back(w1);
                    output.push_back(w2);
                }
                return;
            }
            wrapper_list ch;
            for(int j = 0; j < max_record; j ++)
                ch.push_back(sorted.at(i + j)->assoc_wrapper);
            auto* w = new wrapper;
            set_as_children(w, ch);
            set_bound_rect(w);
            output.push_back(w);
        }
    }
};

template<int _max_record, int _min_record, class _tree>
class rtree_alg_base
{
public:
    typedef _tree mytree;
    typedef typename _tree::value value;
    typedef typename _tree::wrapper wrapper;
    typedef typename _tree::iterator iterator;
    typedef typename _tree::const_iterator const_iterator;
    typedef typename _tree::value value;
    typedef typename value::bind_arg bind_arg;
    typedef rtree_finder<mytree> finder;
    typedef list<value> value_list;
    static const int max_record = _max_record;
    static const int min_record = _min_record;

public:
    rtree_alg_base(mytree& tr): _mytree(tr) {}

protected:
    mytree&         _mytree;

protected:
    static void self_transfer(mytree& mt, iterator i, iterator dp)
    {
        assert(i.is_valid() && dp.is_valid());
        mytree t;
        mt.detach(t, i);
        mt.attach(t, mt.birth_tail(dp));
    }
    static void update_subtree_rect(iterator p)
    {
        assert(p.is_valid() && !p.is_leaf());
        auto i = p.child();
        rectf rc = i->const_rect();
        for(i.to_next(); i.is_valid(); i.to_next())
            union_rect(rc, rc, i->const_rect());
        p->set_rect(rc);
    }
    static void update_rect_recursively(iterator p)
    {
        assert(p.is_valid() && !p.is_leaf());
        auto i = p.child();
        rectf rc = i->const_rect();
        for(i.to_next(); i.is_valid(); i.to_next())
            union_rect(rc, rc, i->const_rect());
        if(p->const_rect() != rc) {
            p->set_rect(rc);
            auto g = p.parent();
            if(g.is_valid())
                update_rect_recursively(g);
        }
    }
    static void update_rect_recursively(iterator p, iterator i)
    {
        assert(p.is_valid() && i.is_valid() && (i.parent() == p));
        rectf rc = p->const_rect();
        union_rect(rc, rc, i->const_rect());
        if(p->const_rect() != rc) {
            p->set_rect(rc);
            auto g = p.parent();
            if(g.is_valid())
                update_rect_recursively(g, p);
        }
    }
    static void collect_leaves(value_list& candidates, iterator i)
    {
        assert(i.is_valid());
        if(i.is_leaf()) {
            candidates.push_back(*i);
            return;
        }
        for(auto j = i.child(); j.is_valid(); j.to_next())
            collect_leaves(candidates, j);
    }

protected:
    void split_binary(iterator p, iterator i, iterator j)
    {
        assert(p.is_valid() && (p.childs() == max_record + 1));
        assert(i.is_valid() && j.is_valid() && (i.parent() == p) && (j.parent() == p));
        auto g = p.parent();
        iterator q;
        if(g.is_valid())
            q = _mytree.insert_after(p);
        else {
            mytree mt;
            auto r = mt.insert(iterator(nullptr));
            r->set_rect(p->const_rect());
            auto p1 = mt.birth(r);
            auto p2 = mt.birth_tail(r);
            mt.attach(_mytree, p1);
            mt.swap(_mytree);
            p = p1;
            q = p2;
        }
        int max_split = max_record + 1 - min_record;
        /* retrieve the split line */
        auto c1 = i->const_rect().center();
        auto c2 = j->const_rect().center();
        float coef[3];
        calc_perpendicular_bisector(coef, c1, c2);
        /* sort them by signed distance */
        typedef multimap<float, iterator> sort_by_sd;
        sort_by_sd sbsd;
        float muler = 1.f / sqrtf(coef[0] * coef[0] + coef[1] * coef[1]);
        for(auto m = p.child(); m.is_valid(); m.to_next()) {
            if(m != i && m != j) {
                float d = muler * calc_orientation(coef, m->const_rect().center());
                sbsd.insert(std::make_pair(d, m));
            }
        }
        /* dump to a random access container */
        typedef vector<iterator> sd_queue;
        assert(sbsd.size() == max_record - 1);
        int divpos = -1, idx = 0;
        sd_queue sdq;
        sdq.resize(sbsd.size());
        for(auto n = sbsd.begin(); n != sbsd.end(); ++ n, ++ idx) {
            sdq.at(idx) = n->second;
            if((divpos < 0) && (n->first >= 0.f))
                divpos = idx;
        }
        if(divpos < 0)
            divpos = 0;
        /* clamp */
        if(divpos < min_record)
            divpos = min_record;
        else if(divpos > max_split)
            divpos = max_split;
        /* transfer from p to q */
        self_transfer(_mytree, j, q);
        for(auto pos = divpos; pos != (int)sdq.size(); ++ pos)
            self_transfer(_mytree, sdq.at(pos), q);
        /* update p, q's rect */
        update_subtree_rect(p);
        update_subtree_rect(q);
    }
};

/*
 * It's a bit weird that the Guttman's paper said that:
 * "PickNext simply chooses any of the remaining entries"
 * Does this algorithm actually useful, or maybe I implemented it wrong.
 * Use quadratic split in preference.
 * plus TODO:
 * R* algorithm to be implemented, but undone.
 */
template<int _max_record, int _min_record, class _tree>
class linear_split_alg:
    public rtree_alg_base<_max_record, _min_record, _tree>
{
public:
    linear_split_alg(mytree& tr):
        rtree_alg_base<_max_record, _min_record, _tree>(tr)
    {
        /* better be static_assert! */
        assert((min_record * 2 <= max_record) &&
            "min_record should be less than half max_record."
            );
    }
    iterator insert(bind_arg ba, const rectf& rc)
    {
        if(!_mytree.is_valid()) {
            auto i = _mytree.insert(iterator(nullptr));
            *i = value(ba, rc);
            return i;
        }
        auto i = find_insert_position(_mytree.get_root());
        if(i.is_valid()) {
            assert(i.childs() < max_record);
            auto j = _mytree.birth_tail(i);
            *j = value(ba, rc);
            update_rect_recursively(i, j);
            return j;
        }
        i = find_first_leaf(_mytree.get_root());
        assert(i.is_valid() && i.is_leaf());
        auto p = i.parent();
        if(!p.is_valid()) {
            assert(i == _mytree.get_root());
            auto j = _mytree.birth(i);
            auto k = _mytree.birth_tail(i);
            *j = *i;
            *k = value(ba, rc);
            rectf new_rc;
            union_rect(new_rc, i->const_rect(), rc);
            i->set_bind_arg(0);
            i->set_rect(new_rc);
            return k;
        }
        assert(p.childs() == max_record);
        /* temporarily insert into p, then split */
        auto j = _mytree.birth_tail(p);
        *j = value(ba, rc);
        update_rect_recursively(p, j);
        do { p = adjust_tree(p); }
        while(p.is_valid());
        return j;
    }
    void remove(bind_arg ba, const rectf& rc)
    {
        auto f = finder(_mytree).find(ba, rc);
        if(!f.is_valid())
            return;
        assert(f.is_leaf());
        if(f.is_root()) {
            _mytree.clear();
            return;
        }
        auto p = f.parent();
        assert(p.is_valid());
        _mytree.erase(f);
        if(p.childs() >= min_record) {
            update_rect_recursively(p);
            return;
        }
        value_list candidates;
        do { p = condense_tree(candidates, p); }
        while(p.is_valid());
        reinsert(candidates);
    }

protected:
    iterator find_insert_position(iterator i)
    {
        assert(i.is_valid());
        if(i.is_leaf())
            return iterator(nullptr);
        auto j = i.child();
        if(j.is_leaf())
            return (i.childs() < max_record) ? i : iterator(nullptr);
        for(; j.is_valid(); j.to_next()) {
            auto r = find_insert_position(j);
            if(r.is_valid())
                return r;
        }
        return iterator(nullptr);
    }
    iterator find_first_leaf(iterator i)
    {
        assert(i.is_valid());
        return i.is_leaf() ? i : find_first_leaf(i.child());
    }
    iterator find_extreme_u(iterator p)
    {
        assert(p.is_valid() && !p.is_leaf());
        auto i = p.child();
        assert(i.is_valid());
        float m = i->const_rect().width();
        auto r = i;
        for(i.to_next(); i.is_valid(); i.to_next()) {
            float u = i->const_rect().width();
            if(u > m)
                r = i, m = u;
        }
        return r;
    }
    iterator find_extreme_v_but_not(iterator p, iterator c)
    {
        assert(p.is_valid() && !p.is_leaf());
        assert(c.is_valid() && (c.parent() == p));
        iterator r;
        float m = -FLT_MAX;
        for(auto i = p.child(); i.is_valid(); i.to_next()) {
            if(i == c)
                continue;
            float v = i->const_rect().height();
            if(v > m)
                r = i, m = v;
        }
        assert(r.is_valid());
        return r;
    }
    void split(iterator p)
    {
        assert(p.is_valid() && (p.childs() == max_record + 1));
        auto i = find_extreme_u(p);
        auto j = find_extreme_v_but_not(p, i);
        split_binary(p, i, j);
    }
    iterator adjust_tree(iterator i)
    {
        assert(i.is_valid() && !i.is_leaf());
        assert(i.childs() == max_record + 1);
        auto p = i.parent();
        split(i);
        if(!p.is_valid())
            return iterator(nullptr);
        return (p.childs() <= max_record) ? iterator(nullptr) : p;
    }
    iterator condense_tree(value_list& candidates, iterator i)
    {
        assert(i.is_valid() && !i.is_leaf());
        assert(i.childs() < min_record);
        if(i.is_root()) {
            if(i.childs() == 1) {
                /* height decrease */
                mytree t;
                _mytree.detach(t, i.child());
                _mytree.swap(t);
            }
            return iterator(nullptr);
        }
        collect_leaves(candidates, i);
        auto p = i.parent();
        assert(p.is_valid());
        _mytree.erase(i);
        if(p.childs() >= min_record) {
            update_subtree_rect(p);
            return iterator(nullptr);
        }
        return p;
    }
    void reinsert(value_list& candidates)
    {
        if(candidates.empty())
            return;
        std::for_each(candidates.begin(), candidates.end(), [this](const value& val) {
            this->insert(val.get_bind_arg(), val.const_rect());
        });
    }
};

template<int _max_record, int _min_record, class _tree>
class quadratic_split_alg:
    public rtree_alg_base<_max_record, _min_record, _tree>
{
public:
    quadratic_split_alg(mytree& tr):
        rtree_alg_base<_max_record, _min_record, _tree>(tr)
    {
        /* better be static_assert! */
        assert((min_record * 2 <= max_record) &&
            "min_record should be less than half max_record."
            );
    }
    iterator insert(bind_arg ba, const rectf& rc)
    {
        if(!_mytree.is_valid()) {
            auto i = _mytree.insert(iterator(nullptr));
            *i = value(ba, rc);
            return i;
        }
        auto i = find_insert_position(_mytree.get_root(), rc);
        assert(i.is_valid());
        auto p = i.parent();
        if(!p.is_valid()) {
            assert(i == _mytree.get_root());
            auto j = _mytree.birth(i);
            auto k = _mytree.birth_tail(i);
            *j = *i;
            *k = value(ba, rc);
            rectf new_rc;
            union_rect(new_rc, i->const_rect(), rc);
            i->set_bind_arg(0);
            i->set_rect(new_rc);
            return k;
        }
        if(p.childs() < max_record) {
            auto j = _mytree.birth_tail(p);
            *j = value(ba, rc);
            update_rect_recursively(p, j);
            return j;
        }
        /* temporarily insert into p, then split */
        auto j = _mytree.birth_tail(p);
        *j = value(ba, rc);
        update_rect_recursively(p, j);
        do { p = adjust_tree(p); }
        while(p.is_valid());
        return j;
    }
    void remove(bind_arg ba, const rectf& rc)
    {
        auto f = finder(_mytree).find(ba, rc);
        if(!f.is_valid())
            return;
        assert(f.is_leaf());
        if(f.is_root()) {
            _mytree.clear();
            return;
        }
        auto p = f.parent();
        assert(p.is_valid());
        _mytree.erase(f);
        if(p.childs() >= min_record) {
            update_rect_recursively(p);
            return;
        }
        value_list candidates;
        do { p = condense_tree(candidates, p); }
        while(p.is_valid());
        reinsert(candidates);
    }

protected:
    iterator find_least_enlargement(iterator i, const rectf& rc)
    {
        assert(i.is_valid());
        auto c = i.child();
        auto ret = c;
        float enlargement = FLT_MAX;
        for(; c.is_valid(); c = c.next()) {
            const rectf& oldrc = c->const_rect();
            float oldarea = oldrc.area();
            rectf fr;
            union_rect(fr, oldrc, rc);
            float area = fr.area();
            float enlarge = area - oldarea;
            if(enlarge < enlargement) {
                enlargement = enlarge;
                ret = c;
            }
        }
        return ret;
    }
    float find_biggest_enlargement(iterator& i, iterator& j, iterator p)
    {
        assert(p.is_valid() && !p.is_leaf());
        assert(p.childs() == max_record + 1);
        iterator f = p.child(), g;
        assert(f.next().is_valid());
        float m = find_biggest_enlargement(g, f, f.next(), p.last_child());
        i = f, j = g;
        for(f.to_next(); f != p.last_child(); f.to_next()) {
            float n = find_biggest_enlargement(g, f, f.next(), p.last_child());
            if(n > m) {
                m = n;
                i = f, j = g;
            }
        }
        return m;
    }
    float find_biggest_enlargement(iterator& j, iterator i, iterator from, iterator to)
    {
        assert(i.is_valid() && from.is_valid() && to.is_valid());
        iterator t = from;
        float s = i->const_rect().area();
        float m = calc_enlarge_area(i->const_rect(), t->const_rect(), s);
        j = t;
        if(t == to)
            return m;
        for(t.to_next();; t.to_next()) {
            float n = calc_enlarge_area(i->const_rect(), t->const_rect(), s);
            if(n > m) { m = n, j = t; }
            if(t == to)
                break;
        }
        return m;
    }
    iterator find_insert_position(iterator i, const rectf& rc)
    {
        if(i.is_leaf())
            return i;
        auto trygetch = i.child();
        if(trygetch.is_leaf())
            return trygetch;
        auto j = find_least_enlargement(i, rc);
        assert(j.is_valid());
        return find_insert_position(j, rc);
    }
    void split(iterator p)
    {
        assert(p.is_valid() && (p.childs() == max_record + 1));
        iterator i, j;
        find_biggest_enlargement(i, j, p);
        split_binary(p, i, j);
    }
    iterator adjust_tree(iterator i)
    {
        assert(i.is_valid() && !i.is_leaf());
        assert(i.childs() == max_record + 1);
        auto p = i.parent();
        split(i);
        if(!p.is_valid())
            return iterator(nullptr);
        return (p.childs() <= max_record) ? iterator(nullptr) : p;
    }
    iterator condense_tree(value_list& candidates, iterator i)
    {
        assert(i.is_valid() && !i.is_leaf());
        assert(i.childs() < min_record);
        if(i.is_root()) {
            if(i.childs() == 1) {
                /* height decrease */
                mytree t;
                _mytree.detach(t, i.child());
                _mytree.swap(t);
            }
            return iterator(nullptr);
        }
        collect_leaves(candidates, i);
        auto p = i.parent();
        assert(p.is_valid());
        _mytree.erase(i);
        if(p.childs() >= min_record) {
            update_subtree_rect(p);
            return iterator(nullptr);
        }
        return p;
    }
    void reinsert(value_list& candidates)
    {
        if(candidates.empty())
            return;
        std::for_each(candidates.begin(), candidates.end(), [this](const value& val) {
            this->insert(val.get_bind_arg(), val.const_rect());
        });
    }
};

template<class _ty,
    class _alg,
    class _wrapper = rtree_node<_ty>,
    class _alloc = _tree_allocator<_wrapper>
    >
class rtree:
    public tree<_ty, _wrapper, _alloc>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _alloc alloc;
    typedef _alg alg;
    typedef rtree<_ty, _alg, _wrapper, _alloc> myref;
    typedef tree<_ty, _wrapper, _alloc> superref;
    typedef typename superref::iterator iterator;
    typedef typename superref::const_iterator const_iterator;
    typedef typename value::bind_arg bind_arg;
    typedef unordered_set<bind_arg> lookup_table;

public:
    template<class _querier, class _cont>
    int query(const _querier& q, _cont& out) const
    {
        if(!is_valid())
            return 0;
        lookup_table lookups;
        auto i = get_root();
        return query_overlapped(i, q, lookups, out);
    }
    template<class _cont>
    int query(const pointf& p1, const pointf& p2, _cont& out) const
    {
        if(!is_valid())
            return 0;
        lookup_table lookups;
        rectf rc;
        rc.set_by_pts(p1, p2);
        auto i = get_root();
        return query_overlapped(i, p1, p2, rc, lookups, out);
    }
    iterator insert(bind_arg ba, const rectf& rc) { return alg(*this).insert(ba, rc); }
    void remove(bind_arg ba, const rectf& rc) { alg(*this).remove(ba, rc); }
    bool empty() const { return !is_valid(); }
    iterator insert_line(bind_arg ba, const pointf& p1, const pointf& p2)
    {
        rectf rc;
        rc.set_by_pts(p1, p2);
        iterator i = insert(ba, rc);
        assert(i);
        i->set_spec(new rtree_spec_line(p1, p2));
        return i;
    }
    template<class _vsl>
    void insert(const _vsl& input)
    {
        if(input.size() == 1) {
            auto& p = input.front();
            insert(p.get_bind_arg(), p.const_rect());
            return;
        }
        rtree_str_load<alg::max_record, alg::min_record, alg::mytree>(*this).load(input);
    }
    void tracing() const
    {
        if(!is_valid())
            return;
        auto r = get_root();
        trace_subtree(r);
        trace_layer(r);
    }

protected:
    template<class _cont>
    int query_overlapped(const_iterator i, const rectf& q, lookup_table& lookups, _cont& out) const
    {
        assert(i.is_valid());
        const auto& rc = i->const_rect();
        if(!is_rect_intersected(rc, q))
            return 0;
        if(i.is_leaf()) {
            if(lookups.find(i->get_bind_arg()) != lookups.end())
                return 0;
            lookups.insert(i->get_bind_arg());
            out.push_back(i->get_bind_arg());
            return 1;
        }
        int c = 0;
        for(auto j = i.child(); j.is_valid(); j = j.next())
            c += query_overlapped(j, q, lookups, out);
        return c;
    }
    template<class _cont>
    int query_overlapped(const_iterator i, const pointf& q, lookup_table& lookups, _cont& out) const
    {
        assert(i.is_valid());
        const auto& rc = i->const_rect();
        if(!rc.in_rect(q))
            return 0;
        if(i.is_leaf()) {
            if(lookups.find(i->get_bind_arg()) != lookups.end())
                return 0;
            lookups.insert(i->get_bind_arg());
            out.push_back(i->get_bind_arg());
            return 1;
        }
        int c = 0;
        for(auto j = i.child(); j.is_valid(); j = j.next())
            c += query_overlapped(j, q, lookups, out);
        return c;
    }
    template<class _cont>
    int query_overlapped(const_iterator i, const pointf& p1, const pointf& p2, const rectf& rc, lookup_table& lookups, _cont& out) const
    {
        assert(i.is_valid());
        if(i.is_leaf()) {
            auto* spec = i->get_spec();
            bool isovl = spec ? spec->query_overlap(p1, p2) : is_rect_intersected(i->const_rect(), rc);
            if(isovl && lookups.find(i->get_bind_arg()) == lookups.end()) {
                lookups.insert(i->get_bind_arg());
                out.push_back(i->get_bind_arg());
                return 1;
            }
            return 0;
        }
        if(!is_rect_intersected(i->const_rect(), rc))
            return 0;
        int c = 0;
        for(auto j = i.child(); j.is_valid(); j = j.next())
            c += query_overlapped(j, p1, p2, rc, lookups, out);
        return c;
    }
    void trace_layer(const_iterator i) const
    {
        assert(i.is_valid());
        int r = rand() % 256;
        int g = rand() % 256;
        int b = rand() % 256;
        string cr;
        cr.format(_t("rgb(%d,%d,%d)"), r, g, b);
        trace(_t("@!\n"));
        trace(_t("@&strokeColor=%s;\n"), cr.c_str());
        trace(_t("@&withArrow=false;\n"));
        for(auto j = i; j.is_valid(); j = j.next()) {
            const auto& rc = j->const_rect();
            trace(_t("@rect %f, %f, %f, %f;\n"), rc.left, rc.top, rc.right, rc.bottom);
        }
        trace(_t("@@\n"));
    }
    void trace_subtree(const_iterator i) const
    {
        assert(i.is_valid());
        if(i.is_leaf())
            return;
        for(auto j = i.child(); j.is_valid(); j = j.next())
            trace_subtree(j);
        trace_layer(i.child());
    }
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef Single_Pin_Subchannel_h
#define Single_Pin_Subchannel_h

#include <vector>

// Trilinos includes
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"

// SCALE includes
#include "Nemesis/harness/DBC.hh"

// vendored includes
#include <gsl/gsl>

namespace enrico {

//===========================================================================//
/*!
 * \class Single_Pin_Subchannel
 * \brief Solve subchannel flow in single channel.
 *
 * This class implements a two-equation subchannel model involving
 * conservation of energy and axial momentum.  No lateral flow is accounted
 * for.  The model excludes friction momentum losses and any spacer grid
 * effects.  It is useful for giving qualitatively-correct behavior but
 * should not be used for actual analysis.
 */
//===========================================================================//
class Single_Pin_Subchannel {
public:
  //@{
  //! Typedefs
  using RCP_PL = Teuchos::RCP<Teuchos::ParameterList>;
  //@}

  enum Verbosity { NONE, LOW, HIGH };

private:
  // >>> DATA

  // Solve parameters
  double d_tol;
  int d_max_iters;
  Verbosity d_verbosity;

  // Axial grid
  std::vector<double> d_delta_z;

  // Channel geometry
  double d_area;

  // Flow conditions
  double d_mdot;
  double d_T_inlet;
  double d_p_exit;

public:
  // Constructor
  Single_Pin_Subchannel(RCP_PL& parameters, const std::vector<double>& delta_z);

  // Set channel area (cm^2)
  void set_channel_area(double area)
  {
    // Channel area internally is needed in m^2
    Expects(area > 0.0);
    d_area = 1e-4 * area;
  }

  // Set mass flow rate (kg/s)
  void set_mass_flow_rate(double mdot)
  {
    Expects(mdot > 0.0);
    Expects(mdot < 2.0);
    d_mdot = mdot;
  }

  // Set inlet temperature (K)
  void set_inlet_temperature(double T)
  {
    Expects(T > 0);
    Expects(T < 1000);
    d_T_inlet = T;
  }

  // Set exit pressure (Pa)
  void set_exit_pressure(double p)
  {
    Expects(p > 0);
    Expects(p < 2.2e7);
    d_p_exit = p;
  }

  // Solve for single subchannel given power distribution
  void solve(const std::vector<double>& power,
             std::vector<double>& temperature,
             std::vector<double>& density);
};

//---------------------------------------------------------------------------//
} // end namespace enrico

//---------------------------------------------------------------------------//
#endif // Single_Pin_Subchannel_h

//---------------------------------------------------------------------------//
// end of Single_Pin_Subchannel.h
//---------------------------------------------------------------------------//

max_stars_repo_name: sphamil/enrico

max_issues_repo_name: sphamil/enrico

max_forks_repo_name: sphamil/enrico


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_syrk (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = -1.0f;
   float beta = 0.1f;
   float A[] = { 0.412f, -0.229f };
   int lda = 1;
   float C[] = { 0.628f, -0.664f, -0.268f, 0.096f };
   int ldc = 2;
   float C_expected[] = { -0.106944f, 0.027948f, -0.268f, -0.042841f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1566)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = -1.0f;
   float beta = 0.1f;
   float A[] = { 0.101f, -0.653f };
   int lda = 2;
   float C[] = { 0.432f, 0.107f, -0.952f, -0.532f };
   int ldc = 2;
   float C_expected[] = { 0.032999f, 0.107f, -0.029247f, -0.479609f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1567)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = 1.0f;
   float beta = 0.1f;
   float A[] = { 0.79f, 0.595f };
   int lda = 2;
   float C[] = { 0.257f, 0.183f, -0.021f, -0.053f };
   int ldc = 2;
   float C_expected[] = { 0.6498f, 0.48835f, -0.021f, 0.348725f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1568)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = 1.0f;
   float beta = 0.1f;
   float A[] = { -0.181f, -0.654f };
   int lda = 1;
   float C[] = { -0.4f, 0.615f, 0.147f, -0.163f };
   int ldc = 2;
   float C_expected[] = { -0.007239f, 0.615f, 0.133074f, 0.411416f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1569)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = 0.0f;
   float beta = -1.0f;
   float A[] = { -0.191f, 0.584f };
   int lda = 1;
   float C[] = { -0.719f, -0.681f, -0.003f, 0.544f };
   int ldc = 2;
   float C_expected[] = { 0.719f, -0.681f, 0.003f, -0.544f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1570)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = 0.0f;
   float beta = -1.0f;
   float A[] = { 0.788f, 0.041f };
   int lda = 2;
   float C[] = { 0.029f, 0.365f, 0.739f, -0.769f };
   int ldc = 2;
   float C_expected[] = { -0.029f, -0.365f, 0.739f, 0.769f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1571)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = -0.3f;
   float beta = -1.0f;
   float A[] = { 0.733f, 0.678f };
   int lda = 2;
   float C[] = { -0.941f, 0.96f, 0.07f, -0.295f };
   int ldc = 2;
   float C_expected[] = { 0.779813f, 0.96f, -0.219092f, 0.157095f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1572)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = -0.3f;
   float beta = -1.0f;
   float A[] = { -0.87f, 0.675f };
   int lda = 1;
   float C[] = { -0.602f, -0.432f, -0.984f, 0.384f };
   int ldc = 2;
   float C_expected[] = { 0.37493f, 0.608175f, -0.984f, -0.520687f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1573)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 0.1;
   double beta = -0.3;
   double A[] = { 0.169, -0.875 };
   int lda = 1;
   double C[] = { 0.159, 0.277, 0.865, 0.346 };
   int ldc = 2;
   double C_expected[] = { -0.0448439, -0.0978875, 0.865, -0.0272375 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1574)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 0.1;
   double beta = -0.3;
   double A[] = { 0.536, -0.725 };
   int lda = 2;
   double C[] = { 0.154, -0.445, -0.841, -0.91 };
   int ldc = 2;
   double C_expected[] = { -0.0174704, -0.445, 0.21344, 0.3255625 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1575)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = 0;
   double beta = -1;
   double A[] = { -0.07, 0.8 };
   int lda = 2;
   double C[] = { 0.823, -0.88, -0.136, 0.793 };
   int ldc = 2;
   double C_expected[] = { -0.823, 0.88, -0.136, -0.793 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1576)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = 0;
   double beta = -1;
   double A[] = { -0.058, 0.649 };
   int lda = 1;
   double C[] = { -0.187, 0.294, -0.004, -0.933 };
   int ldc = 2;
   double C_expected[] = { 0.187, 0.294, 0.004, 0.933 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1577)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 1;
   double beta = -1;
   double A[] = { 0.263, -0.289 };
   int lda = 1;
   double C[] = { 0.554, -0.679, 0.993, 0.758 };
   int ldc = 2;
   double C_expected[] = { -0.484831, -0.679, -1.069007, -0.674479 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1578)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 1;
   double beta = -1;
   double A[] = { -0.265, -0.837 };
   int lda = 2;
   double C[] = { -0.994, 0.967, -0.34, -0.069 };
   int ldc = 2;
   double C_expected[] = { 1.064225, -0.745195, -0.34, 0.769569 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1579)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { -0.464, 0.394 };
   int lda = 2;
   double C[] = { -0.45, -0.447, 0.649, 0.055 };
   int ldc = 2;
   double C_expected[] = { -0.5145888, -0.447, 0.7038448, 0.0084292 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1580)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { 0.815, 0.168 };
   int lda = 1;
   double C[] = { 0.817, -0.957, -0.395, -0.382 };
   int ldc = 2;
   double C_expected[] = { 0.6177325, -0.998076, -0.395, -0.3904672 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1581)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.447f, -0.507f, -0.425f, 0.701f };
   int lda = 1;
   float C[] = { 0.16f, -0.245f, 0.922f, -0.437f, 0.24f, 0.008f, -0.095f, 0.749f };
   int ldc = 2;
   float C_expected[] = { -0.0235f, 0.0895f, -0.2329f, 0.2233f, 0.24f, 0.008f, -0.0464f, -0.2342f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1582) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1582) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { -0.421f, -0.435f, -0.914f, -0.493f };
   int lda = 2;
   float C[] = { -0.761f, -0.38f, 0.043f, -0.999f, 0.779f, 0.238f, 0.082f, 0.394f };
   int ldc = 2;
   float C_expected[] = { 0.2663f, 0.0379f, 0.043f, -0.999f, -0.2575f, 0.0065f, -0.064f, -0.11f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1583) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1583) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.827f, -0.896f, 0.417f, 0.865f };
   int lda = 2;
   float C[] = { -0.349f, -0.31f, 0.972f, 0.794f, -0.906f, -0.595f, -0.089f, -0.333f };
   int ldc = 2;
   float C_expected[] = { 0.254587f, 1.54008f, -1.4909f, -0.482723f, -0.906f, -0.595f, 0.634336f, -0.63041f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1584) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1584) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.607f, 0.747f, -0.889f, 0.333f };
   int lda = 1;
   float C[] = { 0.244f, 0.564f, 0.009f, 0.578f, -0.827f, 0.558f, -0.337f, 0.731f };
   int ldc = 2;
   float C_expected[] = { 0.05996f, -1.05166f, 0.009f, 0.578f, 0.980674f, 0.211852f, -0.651432f, 0.339074f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1585) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1585) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { 0.784f, -0.281f, -0.88f, 0.479f };
   int lda = 1;
   float C[] = { 0.491f, 0.531f, 0.805f, -0.097f, 0.728f, 0.674f, -0.705f, -0.754f };
   int ldc = 2;
   float C_expected[] = { 0.004695f, 0.050392f, 0.805f, -0.097f, -1.22932f, 1.35082f, 1.29896f, -1.54804f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1586) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1586) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { 0.272f, -0.146f, 0.155f, 0.038f };
   int lda = 2;
   float C[] = { 0.533f, -0.41f, -0.904f, 0.301f, -0.836f, 0.57f, -0.374f, -0.293f };
   int ldc = 2;
   float C_expected[] = { 0.462668f, 0.453576f, -0.253292f, -0.916294f, -0.836f, 0.57f, 0.315581f, -0.36222f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1587) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1587) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { -0.055f, -0.127f, -0.896f, -0.625f };
   int lda = 2;
   float C[] = { -0.619f, 0.511f, -0.877f, 0.557f, -0.801f, -0.437f, -0.922f, 0.332f };
   int ldc = 2;
   float C_expected[] = { 0.60503f, -0.524104f, -0.877f, 0.557f, 0.652833f, 0.406905f, -0.198f, 0.080191f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1588) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1588) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { -0.528f, 0.759f, -0.079f, 0.952f };
   int lda = 1;
   float C[] = { 0.775f, 0.855f, 0.786f, 0.525f, 0.85f, 0.044f, 0.658f, 0.947f };
   int ldc = 2;
   float C_expected[] = { 0.026504f, -1.1523f, -0.223383f, -1.20586f, 0.85f, 0.044f, -0.507584f, -1.84706f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1589) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1589) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {1, 0};
   double beta[2] = {1, 0};
   double A[] = { -0.049, -0.687, -0.434, 0.294 };
   int lda = 1;
   double C[] = { 0.937, -0.113, 0.796, 0.293, 0.876, -0.199, -0.757, -0.103 };
   int ldc = 2;
   double C_expected[] = { 0.467432, -0.045674, 1.019244, 0.576752, 0.876, -0.199, -0.65508, -0.358192 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1590) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1590) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {1, 0};
   double beta[2] = {1, 0};
   double A[] = { 0.359, -0.364, 0.926, -0.69 };
   int lda = 2;
   double C[] = { 0.306, 0.249, 0.28, 0.229, 0.866, 0.092, 0.886, -0.283 };
   int ldc = 2;
   double C_expected[] = { 0.302385, -0.012352, 0.28, 0.229, 0.947274, -0.492774, 1.267376, -1.56088 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1591) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1591) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {0, 0};
   double A[] = { 0.607, 0.555, -0.85, 0.831 };
   int lda = 2;
   double C[] = { 0.069, 0.368, 0.551, -0.912, -0.243, -0.063, -0.924, 0.192 };
   int ldc = 2;
   double C_expected[] = { -0.0855042, -0.1960886, 0.2898798, -0.1075156, -0.243, -0.063, 0.1316883, 0.4270039 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1592) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1592) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {0, 0};
   double A[] = { 0.427, 0.86, -0.136, 0.002 };
   int lda = 1;
   double C[] = { 0.398, -0.47, 0.011, -0.547, -0.106, 0.016, 0.681, 0.246 };
   int ldc = 2;
   double C_expected[] = { 0.0937373, -0.2760591, 0.011, -0.547, 0.0295482, 0.0288526, -0.0054932, 0.0020124 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1593) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1593) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { 0.718, 0.023, 0.355, -0.492 };
   int lda = 1;
   double C[] = { -0.637, -0.727, -0.475, -0.776, 0.802, -0.55, -0.837, 0.222 };
   int ldc = 2;
   double C_expected[] = { -0.7948013, -0.6854089, -0.475, -0.776, 0.7566473, -0.4198521, -0.7672563, 0.3151921 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1594) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1594) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { 0.209, 0.139, -0.202, -0.223 };
   int lda = 2;
   double C[] = { -0.695, 0.524, 0.212, -0.88, -0.752, 0.291, 0.684, -0.124 };
   int ldc = 2;
   double C_expected[] = { -0.7081182, 0.5090054, 0.2228348, -0.8587166, -0.752, 0.291, 0.6776683, -0.1519201 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1595) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1595) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.365, -0.624, 0.632, 0.348 };
   int lda = 2;
   double C[] = { 0.877, 0.927, -0.377, 0.967, 0.008, 0.292, -0.779, 0.794 };
   int ldc = 2;
   double C_expected[] = { 0.9082933, 0.7647289, -0.377, 0.967, 0.0641972, 0.4470636, -0.9064832, 0.6898704 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1596) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1596) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.067, -0.586, 0.208, 0.331 };
   int lda = 1;
   double C[] = { 0.584, -0.454, 0.93, 0.782, 0.489, -0.278, 0.081, -0.919 };
   int ldc = 2;
   double C_expected[] = { 0.6778197, -0.5114479, 0.8903975, 0.8432225, 0.489, -0.278, 0.0871195, -0.9669385 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1597) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1597) imag");
     };
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc