text:
//
//  her.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 12/23/11.
//  Copyright (c) 2011 University of Colorado Denver. All rights reserved.
//

#ifndef _her_h
#define _her_h

/// @file her.h Performs complex vector outer product.

#include <cctype>
#include "latl.h"

namespace LATL
{
   /// @brief Performs a vector outer product of a complex vector with itself.
   ///
   /// For a complex Hermitian matrix A, complex vector x and real scalar alpha,
   ///
   ///        A := alpha * x * x' + A
   ///
   /// is computed.  The result is Hermitian and is stored as either upper or lower triangular in A.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param uplo Specifies whether A is stored as upper or lower triangular.
   ///
   ///        if uplo = 'U' or 'u' then A is upper triangular
   ///        if uplo = 'L' or 'l' then A is lower triangular
   ///
   /// @param n Specifies the order of the matrix A.  n>=0
   /// @param alpha Real scalar.
   /// @param x Pointer to complex vector x.
   /// @param incx Increment of the vector x.  x!=0
   /// @param A Pointer to complex Hermitian n-by-n matrix A.  If uplo = 'U' or 'u' then only the upper
   /// triangular part of A is referenced and the lower part is not referenced.  If uplo = 'L' or 'l' then
   /// only the lower triangular part of A is referenced and the upper part is not referenced.
   /// @param ldA Column length of matrix A.  ldA>=n.
   /// @ingroup BLAS

   template <typename real_t>
   int HER(char uplo, int_t n, real_t alpha, complex<real_t> *x, int_t incx, complex<real_t> *A, int_t ldA)
   {
      using std::conj;
      using std::real;
      using std::toupper;
      const real_t zero(0.0);
      complex<real_t> t;
      int_t i,j,kx,jx,ix;

      uplo=toupper(uplo);

      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(ldA<n)
         return -7;
      else if((n==0)||(alpha==zero))
         return 0;

      if(incx==1)
      {
         if(uplo=='U')
         {
            for(j=0;j<n;j++)
            {
               t=alpha*conj(x[j]);
               for(i=0;i<j;i++)
                  A[i]+=x[i]*t;
               A[j]=real(A[j])+real(x[j]*t);
               A+=ldA;
            }
         }
         else
         {
            for(j=0;j<n;j++)
            {
               t=alpha*conj(x[j]);
               A[j]=real(A[j])+real(t*x[j]);
               for(i=j+1;i<n;i++)
                  A[i]+=x[i]*t;
               A+=ldA;
            }
         }
      }
      else
      {
         kx=(incx>0)?0:(1-n)*incx;
         if(uplo=='U')
         {
            jx=kx;
            for(j=0;j<n;j++)
            {
               t=alpha*conj(x[jx]);
               ix=kx;
               for(i=0;i<j;i++)
               {
                  A[i]+=x[ix]*t;
                  ix+=incx;
               }
               A[j]=real(A[j])+real(x[jx]*t);
               A+=ldA;
               jx+=incx;
            }
         }
         else
         {
            jx=kx;
            for(j=0;j<n;j++)
            {
               t=alpha*conj(x[jx]);
               A[j]=real(A[j])+real(t*x[jx]);
               ix=jx;
               for(i=j+1;i<n;i++)
               {
                  ix+=incx;
                  A[i]+=x[ix]*t;
               }
               jx+=incx;
               A+=ldA;
            }
         }
      }
      return 0;
   }

#ifdef __latl_cblas
#include <cblas.h>

   template <> int HER<float>(char uplo, int_t n, float alpha, complex<float> *x, int_t incx, complex<float> *A, int_t ldA)
   {
      uplo=std::toupper(uplo);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(ldA<n)
         return -7;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_cher(CblasColMajor,Uplo,n,alpha,x,incx,A,ldA);

      return 0;
   }

   template <> int HER<double>(char uplo, int_t n, double alpha, complex<double> *x, int_t incx, complex<double> *A, int_t ldA)
   {
      uplo=std::toupper(uplo);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(ldA<n)
         return -7;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_zher(CblasColMajor,Uplo,n,alpha,x,incx,A,ldA);

      return 0;
   }
#endif
}
#endif

ext: h

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
// many-small-dgemms.c:
// Mini-app which investigates different ways to perform many small dgemm operations.
// Looping over thousands of rocblas calls is slow, so looking into other solutions.
// This functionality is very interesting to several applications, including nuccor.
// Date: Oct 15th, 2019
// Author: Justin Gage Lietz
// Contact: lietzjg@ornl.gov

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cblas.h>
/* #include <cublas_v2.h> */
#include <hip/hip_runtime.h>
#include <rocblas.h>
#include <omp.h>
#include <string.h> //for memcpy

#include <magma_v2.h>

// tolerance as 1.e-15 fails
#define TOL 1.e-14

// Pulled from magma test code
#define TESTING_CHECK( err )                                                 \
    do {                                                                     \
        magma_int_t err_ = (err);                                            \
        if ( err_ != 0 ) {                                                   \
            fprintf( stderr, "Error: %s\nfailed at %s:%d: error %lld: %s\n", \
                     #err, __FILE__, __LINE__,                               \
                     (long long) err_, magma_strerror(err_) );               \
            exit(1);                                                         \
        }                                                                    \
    } while( 0 )

size_t idx(size_t i, size_t j, size_t nrows, size_t ncols){
  return j*nrows + i;
}

size_t cidx(size_t i, size_t j, size_t nrows, size_t ncols){
  return i*ncols + j;
}

void print_c_matrix(double *mat, size_t nrows, size_t ncols){
  for(int i = 0; i < nrows; i++){
    for(int j = 0; j < ncols; j++){
      printf("%f ", mat[ cidx(i,j,nrows,ncols) ]);
    }
    printf("\n");
  }
  printf("\n");
}

void print_matrix(double *mat, size_t nrows, size_t ncols){
  for(int i = 0; i < nrows; i++){
    for(int j = 0; j < ncols; j++){
      printf("%f ", mat[ idx(i,j,nrows,ncols) ]);
    }
    printf("\n");
  }
  printf("\n");
}

void print_matrix_array(double *mat_A, size_t *Aoffsets, size_t *nrows, size_t *ncols, size_t nBlocks){
  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    size_t Aoffset = Aoffsets[iBlock];
    printf("iBlock: %zu\n",iBlock);
    print_matrix( &mat_A[Aoffset], nrows[iBlock], ncols[iBlock] );
  }
}

void cblas_wrapper(double *A, double* B, double* C, size_t M_in, size_t N_in, size_t K_in){
  int M = M_in;
  int K = K_in;
  int N = N_in;
  int LDA = M;
  int LDB = K;
  int LDC = M;
  double alpha = 1.0;
  double beta = 0.0;

  cblas_dgemm(/*CBLAS*/ CblasColMajor,
      /* TRANS A */ CblasNoTrans,
      /* TRANS B */ CblasNoTrans,
      /* M */M,
      /* N */N,
      /* K */K,
      /* alpha */alpha,
      /* A */A,
      /* LDA */LDA,
      /* B  */B,
      /* LDB */LDB,
      /* BETA */beta,
      /* C */C,
      /* LDC */LDC);
}

void device_dgemm(rocblas_handle handle, double *A, double* B, double* C, size_t M_in, size_t N_in, size_t K_in){
  int M = M_in;
  int K = K_in;
  int N = N_in;
  int LDA = M;
  int LDB = K;
  int LDC = M;
  double alpha = 1.0;
  double beta = 0.0;

  size_t ABytes = sizeof(double)*M*K;
  size_t BBytes = sizeof(double)*K*N;
  size_t CBytes = sizeof(double)*M*N;

  double *d_A;
  double *d_B;
  double *d_C;
  const rocblas_operation transA = rocblas_operation_none;
  const rocblas_operation transB = rocblas_operation_none;

  hipMalloc(&d_A, ABytes);
  hipMalloc(&d_B, BBytes);
  hipMalloc(&d_C, CBytes);

  hipMemcpy(d_A, A, ABytes, hipMemcpyHostToDevice);
  hipMemcpy(d_B, B, BBytes, hipMemcpyHostToDevice);

  rocblas_dgemm(/*rocblas handle*/ handle,
      /* TRANS A */ transA,
      /* TRANS B */ transB,
      /* M */M,
      /* N */N,
      /* K */K,
      /* alpha */&alpha,
      /* A */d_A,
      /* LDA */LDA,
      /* B  */d_B,
      /* LDB */LDB,
      /* BETA */&beta,
      /* C */d_C,
      /* LDC */LDC);

  hipMemcpy(C, d_C, CBytes, hipMemcpyDeviceToHost);

  hipFree(d_A);
  hipFree(d_B);
  hipFree(d_C);
}

void rocblas_array_dgemm(
    rocblas_handle handle,
    const double *A, const double* B, double* C,
    const size_t *Ms, const size_t *Ns, const size_t *Ks,
    const size_t *Aoffsets, const size_t *Boffsets, const size_t *Coffsets, const size_t nBlocks){

  const rocblas_operation transA = rocblas_operation_none;
  const rocblas_operation transB = rocblas_operation_none;

  double *d_A, *d_B, *d_C;
  size_t ABytes = Aoffsets[nBlocks]*sizeof(double);
  size_t BBytes = Boffsets[nBlocks]*sizeof(double);
  size_t CBytes = Coffsets[nBlocks]*sizeof(double);

  hipMalloc(&d_A, ABytes);
  hipMalloc(&d_B, BBytes);
  hipMalloc(&d_C, CBytes);

  hipMemcpy(d_A, A, ABytes, hipMemcpyHostToDevice);
  hipMemcpy(d_B, B, BBytes, hipMemcpyHostToDevice);

  double start = omp_get_wtime();

  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    size_t Aoffset = Aoffsets[iBlock];
    size_t Boffset = Boffsets[iBlock];
    size_t Coffset = Coffsets[iBlock];

    int M = Ms[iBlock];
    int N = Ns[iBlock];
    int K = Ks[iBlock];
    int LDA = M;
    int LDB = K;
    int LDC = M;
    double alpha = 1.0;
    double beta = 0.0;
    /* char transa = 'N'; */
    /* char transb = 'N'; */

    rocblas_dgemm(/*rocblas handle*/ handle,
        /* TRANS A */ transA,
        /* TRANS B */ transB,
        /* M */       M,
        /* N */       N,
        /* K */       K,
        /* alpha */   &alpha,
        /* A */       &d_A[Aoffset],
        /* LDA */     LDA,
        /* B  */      &d_B[Boffset],
        /* LDB */     LDB,
        /* BETA */    &beta,
        /* C */       &d_C[Coffset],
        /* LDC */     LDC);
  }
  double stop = omp_get_wtime();
  printf("time of rocblas_dgemm loop: %f\n", stop - start);
  hipMemcpy(C, d_C, CBytes, hipMemcpyDeviceToHost);

  hipFree(d_A);
  hipFree(d_B);
  hipFree(d_C);
}

void magma_array_dgemm(
    const double *A, const double* B, double* C,
    const size_t *Ms, const size_t *Ns, const size_t *Ks,
    const size_t *Aoffsets, const size_t *Boffsets, const size_t *Coffsets, const size_t nBlocks){
  magma_trans_t transA = MagmaNoTrans;
  magma_trans_t transB = MagmaNoTrans;
  double const* * dA_array;
  double const* * dB_array;
  double **dC_array;
  double 	alpha;
  double 	beta;
  magma_int_t* d_lddb;
  magma_int_t* d_ldda;
  magma_int_t* d_lddc;
  magma_int_t* d_m;
  magma_int_t* d_n;
  magma_int_t* d_k;
  int* h_m;
  int* h_n;
  int* h_k;

  double **hA_array;
  double **hB_array;
  double **hC_array;
  double *d_A_elems;
  double *d_B_elems;
  double *d_C_elems;

  magma_int_t batchCount;
  magma_queue_t queue;
  magma_device_t device;

  magma_getdevice( &device );
  magma_queue_create( device, &queue );

  batchCount = nBlocks;
  // Magma needs larger arrays
  TESTING_CHECK( magma_malloc((void**)&d_m, (batchCount+1)*sizeof(magma_int_t)) );
  TESTING_CHECK( magma_malloc((void**)&d_n, (batchCount+1)*sizeof(magma_int_t)) );
  TESTING_CHECK( magma_malloc((void**)&d_k, (batchCount+1)*sizeof(magma_int_t)) );

  TESTING_CHECK( magma_malloc((void**)&d_ldda, (batchCount+1)*sizeof(magma_int_t) ) );
  TESTING_CHECK( magma_malloc((void**)&d_lddb, (batchCount+1)*sizeof(magma_int_t) ) );
  TESTING_CHECK( magma_malloc((void**)&d_lddc, (batchCount+1)*sizeof(magma_int_t) ) );

  // Looks like magma_malloc works with regular c int as well.
  TESTING_CHECK( magma_malloc_cpu((void**)&h_m, (batchCount+1)*sizeof(int) ) );
  TESTING_CHECK( magma_malloc_cpu((void**)&h_n, (batchCount+1)*sizeof(int) ) );
  TESTING_CHECK( magma_malloc_cpu((void**)&h_k, (batchCount+1)*sizeof(int) ) );

  // dA_array is the array of pointers need by dgemm
  // d_A_elems are the actual mtx elements being pointed to
  // hA_array is the host side pointers that will get passed to dA_array
  TESTING_CHECK( magma_malloc( (void**)&dA_array, sizeof(double*)*batchCount ) );
  TESTING_CHECK( magma_malloc( (void**)&dB_array, sizeof(double*)*batchCount ) );
  TESTING_CHECK( magma_malloc( (void**)&dC_array, sizeof(double*)*batchCount ) );

  TESTING_CHECK( magma_malloc( (void**)&d_A_elems, sizeof(double)*Aoffsets[nBlocks] ) );
  TESTING_CHECK( magma_malloc( (void**)&d_B_elems, sizeof(double)*Boffsets[nBlocks] ) );
  TESTING_CHECK( magma_malloc( (void**)&d_C_elems, sizeof(double)*Coffsets[nBlocks] ) );

  TESTING_CHECK( magma_malloc_cpu( (void**)&hA_array, sizeof(double*)*batchCount ) );
  TESTING_CHECK( magma_malloc_cpu( (void**)&hB_array, sizeof(double*)*batchCount ) );
  TESTING_CHECK( magma_malloc_cpu( (void**)&hC_array, sizeof(double*)*batchCount ) );

  // magma needs 32-bit int, while main was using 64
  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    h_m[iBlock] = Ms[iBlock];
    h_n[iBlock] = Ns[iBlock];
    h_k[iBlock] = Ks[iBlock];
  }

  magma_setvector(batchCount, sizeof(magma_int_t), h_m, 1, d_m, 1, queue);
  magma_setvector(batchCount, sizeof(magma_int_t), h_n, 1, d_n, 1, queue);
  magma_setvector(batchCount, sizeof(magma_int_t), h_k, 1, d_k, 1, queue);
  magma_setvector(batchCount, sizeof(magma_int_t), h_m, 1, d_ldda, 1, queue);
  magma_setvector(batchCount, sizeof(magma_int_t), h_k, 1, d_lddb, 1, queue);
  magma_setvector(batchCount, sizeof(magma_int_t), h_m, 1, d_lddc, 1, queue);

  magma_setvector(Aoffsets[nBlocks], sizeof(double), A, 1, d_A_elems, 1, queue);
  magma_setvector(Boffsets[nBlocks], sizeof(double), B, 1, d_B_elems, 1, queue);

  size_t Aoffset;
  size_t Boffset;
  size_t Coffset;

  for(int iBlock = 0; iBlock < nBlocks; iBlock++){
    Aoffset = Aoffsets[iBlock];
    Boffset = Boffsets[iBlock];
    Coffset = Coffsets[iBlock];

    hA_array[iBlock] = d_A_elems + Aoffset;
    hB_array[iBlock] = d_B_elems + Boffset;
    hC_array[iBlock] = d_C_elems + Coffset;
  }

  magma_setvector(batchCount, sizeof(double*), hA_array, 1, dA_array, 1, queue);
  magma_setvector(batchCount, sizeof(double*), hB_array, 1, dB_array, 1, queue);
  magma_setvector(batchCount, sizeof(double*), hC_array, 1, dC_array, 1, queue);

  alpha = 1.0;
  beta = 0.0;
  batchCount = nBlocks;

  double start,stop;

  start = omp_get_wtime();

  magmablas_dgemm_vbatched(	      transA,
      /* magma_trans_t */ 	      transB,
      /* magma_int_t * */         d_m,
      /* magma_int_t * */	        d_n,
      /* magma_int_t * */	        d_k,
      /* double */	              alpha,
      /* double const *const * */	dA_array,
      /* magma_int_t * */	        d_ldda,
      /* double const *const * */	dB_array,
      /* magma_int_t * */	        d_lddb,
      /* double */	              beta,
      /* double ** */	            dC_array,
      /* magma_int_t * */	        d_lddc,
      /* magma_int_t */	          batchCount,
      /* magma_queue_t */	        queue);
  stop = omp_get_wtime();

  printf("time of magmablas_dgemm_vbatched: %f\n", stop - start);

  magma_getvector(Coffsets[nBlocks], sizeof(double), d_C_elems, 1, C, 1, queue);

  // Clean up
  magma_queue_destroy( queue );

  TESTING_CHECK( magma_free(d_m) );
  TESTING_CHECK( magma_free(d_n) );
  TESTING_CHECK( magma_free(d_k) );

  TESTING_CHECK( magma_free(d_ldda) );
  TESTING_CHECK( magma_free(d_lddb) );
  TESTING_CHECK( magma_free(d_lddc) );

  TESTING_CHECK( magma_free_cpu(h_m) );
  TESTING_CHECK( magma_free_cpu(h_n) );
  TESTING_CHECK( magma_free_cpu(h_k) );

  TESTING_CHECK( magma_free(dA_array) );
  TESTING_CHECK( magma_free(dB_array) );
  TESTING_CHECK( magma_free(dC_array) );

  TESTING_CHECK( magma_free(d_A_elems) );
  TESTING_CHECK( magma_free(d_B_elems) );
  TESTING_CHECK( magma_free(d_C_elems) );

  TESTING_CHECK( magma_free_cpu(hA_array) );
  TESTING_CHECK( magma_free_cpu(hB_array) );
  TESTING_CHECK( magma_free_cpu(hC_array) );
}

void host_matmul(const double *A, const double* B, double* C, const size_t M, const size_t N, const size_t K){
  for(size_t i = 0; i < M; i++){
    for(size_t j = 0; j < N; j++){
      C[ idx(i,j,M,N) ] = 0.0;
      for(size_t k = 0; k < K; k++){
        C[ idx(i,j,M,N) ] += A[ idx(i,k,M,K) ] * B[ idx(k,j,K,N) ];
      }
    }
  }
}

void host_array_matmul(const double *A, const double* B, double* C,
    const size_t *Ms, const size_t *Ns, const size_t *Ks,
    const size_t *Aoffsets, const size_t *Boffsets, const size_t *Coffsets, const size_t nBlocks){

  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    size_t Aoffset = Aoffsets[iBlock];
    size_t Boffset = Boffsets[iBlock];
    size_t Coffset = Coffsets[iBlock];

    host_matmul(&A[Aoffset], &B[Boffset], &C[Coffset], Ms[iBlock], Ns[iBlock], Ks[iBlock]);
  }
}

void host_array_dgemm(const double *A, const double* B, double* C,
    const size_t *Ms, const size_t *Ns, const size_t *Ks,
    const size_t *Aoffsets, const size_t *Boffsets, const size_t *Coffsets, const size_t nBlocks){

  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    size_t Aoffset = Aoffsets[iBlock];
    size_t Boffset = Boffsets[iBlock];
    size_t Coffset = Coffsets[iBlock];

    int M = Ms[iBlock];
    int N = Ns[iBlock];
    int K = Ks[iBlock];
    int LDA = M;
    int LDB = K;
    int LDC = M;
    double alpha = 1.0;
    double beta = 0.0;
    /* char transa = 'N'; */
    /* char transb = 'N'; */

    cblas_dgemm(/*CBLAS*/ CblasColMajor,
        /* TRANS A */ CblasNoTrans,
        /* TRANS B */ CblasNoTrans,
        /* M */       M,
        /* N */       N,
        /* K */       K,
        /* alpha */   alpha,
        /* A */       &A[Aoffset],
        /* LDA */     LDA,
        /* B  */      &B[Boffset],
        /* LDB */     LDB,
        /* BETA */    beta,
        /* C */       &C[Coffset],
        /* LDC */     LDC);
  }
}

void host_transpose(double *mat_in, double *mat_out, size_t nrows_in, size_t ncols_in){
  for(size_t i = 0; i < nrows_in; i++){
    for(size_t j = 0; j < ncols_in; j++){
      mat_out[ idx(j,i,ncols_in,nrows_in) ] = mat_in[ idx(i,j,nrows_in, ncols_in) ];
    }
  }
}

int check_equal(double *mat_A, double *mat_B, size_t nrows, size_t ncols){
  int flag = 1;
  for(size_t i = 0; i < nrows; i++){
    for(size_t j = 0; j < ncols; j++){
      if( fabs( mat_A[ idx(i,j,nrows,ncols) ] - mat_B[ idx(i,j,nrows,ncols) ]) > TOL ){
        flag = 0;
        return flag;
      }
    }
  }
  return flag;
}

int check_matrix_array_equal(double *mat_A, double *mat_B, size_t *Aoffsets, size_t *Boffsets, size_t *nrows, size_t *ncols, size_t nBlocks){
  int flag = 1;
  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    size_t Aoffset = Aoffsets[iBlock];
    size_t Boffset = Boffsets[iBlock];

    flag *= check_equal( &mat_A[Aoffset], &mat_B[Boffset], nrows[iBlock], ncols[iBlock] );
    if( !check_equal(&mat_A[Aoffset], &mat_B[Boffset], nrows[iBlock], ncols[iBlock]) ){
      printf("error in block: %zu\n");
      printf("A\n");
      print_matrix( &mat_A[Aoffset], nrows[iBlock], ncols[iBlock] );
      printf("B\n");
      print_matrix( &mat_B[Aoffset], nrows[iBlock], ncols[iBlock] );
      return flag;
    }
  }
  return flag;
}

void generate_random_array(size_t* array, size_t nBlocks, size_t max_size){
  size_t rando;

  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    rando = rand()%max_size + 1;
    array[iBlock] = rando;
  }
}

void load_random_matrix(double* mat_in, size_t nrows, size_t ncols){
  double rando;

  for(size_t i = 0; i < nrows; i++){
    for(size_t j = 0; j < ncols; j++){
      rando= (double) rand()/RAND_MAX*2.0 - 1.0;
      mat_in[ idx(i,j,nrows,ncols) ] = rando;
    }
  }
}

int main(){
  // Start magma
  TESTING_CHECK( magma_init() );
  magma_print_environment();

  // magma seems to break at 2^16, but 2^15 blocks is fine
  size_t nBlocks = 1<<15;
  size_t max_block_size = 10;
  // Small test for debugging
  // size_t nBlocks = 3;
  // size_t max_block_size = 4;

  size_t *Ms_array;
  size_t *Ns_array;
  size_t *Ks_array;
  size_t *Aoffsets;
  size_t *Boffsets;
  size_t *Coffsets;

  Ms_array = (size_t*)malloc(sizeof(size_t)*nBlocks);
  Ns_array = (size_t*)malloc(sizeof(size_t)*nBlocks);
  Ks_array = (size_t*)malloc(sizeof(size_t)*nBlocks);

  Aoffsets = (size_t*)malloc(sizeof(size_t)*(nBlocks+1));
  Boffsets = (size_t*)malloc(sizeof(size_t)*(nBlocks+1));
  Coffsets = (size_t*)malloc(sizeof(size_t)*(nBlocks+1));

  /* srand( time(NULL) ); */
  // Reproducible seed
  srand( 100 );

  generate_random_array(Ms_array,nBlocks,max_block_size);
  generate_random_array(Ns_array,nBlocks,max_block_size);
  generate_random_array(Ks_array,nBlocks,max_block_size);

  Aoffsets[0] = 0;
  Boffsets[0] = 0;
  Coffsets[0] = 0;

  // Offset array elements indicate starting location of each block
  // Final array element indicates where the next block WOULD start
  // This is needed to count the total number of elements
  for(size_t iBlock = 1; iBlock < nBlocks + 1; iBlock++){
    Aoffsets[iBlock] = Aoffsets[iBlock-1] + Ms_array[iBlock-1] * Ks_array[iBlock-1];
    Boffsets[iBlock] = Boffsets[iBlock-1] + Ks_array[iBlock-1] * Ns_array[iBlock-1];
    Coffsets[iBlock] = Coffsets[iBlock-1] + Ms_array[iBlock-1] * Ns_array[iBlock-1];
  }

  size_t ABytes = Aoffsets[nBlocks]*sizeof(double);
  size_t BBytes = Boffsets[nBlocks]*sizeof(double);
  size_t CBytes = Coffsets[nBlocks]*sizeof(double);

  double *h_A;
  double *h_B;
  double *h_C;
  double *h_C_host;
  double *h_C_device;
  double *h_C_magma;

  h_A = (double*)malloc(ABytes);
  h_B = (double*)malloc(BBytes);
  h_C = (double*)malloc(CBytes);
  h_C_host = (double*)malloc(CBytes);
  h_C_device = (double*)malloc(CBytes);
  h_C_magma = (double*)malloc(CBytes);

  for(size_t iBlock = 0; iBlock < nBlocks; iBlock++){
    size_t Aoffset = Aoffsets[iBlock];
    size_t Boffset = Boffsets[iBlock];

    load_random_matrix(&h_A[Aoffset], Ms_array[iBlock], Ks_array[iBlock]);
    load_random_matrix(&h_B[Boffset], Ks_array[iBlock], Ns_array[iBlock]);
  }

  double start,stop;
  int check;
  //////////////////////////////////////////////////////////////////////
  // Tests on host
  //////////////////////////////////////////////////////////////////////

  start = omp_get_wtime();
  host_array_matmul(h_A, h_B, h_C, Ms_array, Ns_array, Ks_array, Aoffsets, Boffsets, Coffsets, nBlocks);
  stop = omp_get_wtime();
  printf("naive loop time: %f\n", stop - start);

  start = omp_get_wtime();
  host_array_dgemm(h_A, h_B, h_C_host, Ms_array, Ns_array, Ks_array, Aoffsets, Boffsets, Coffsets, nBlocks);
  stop = omp_get_wtime();
  printf("host blas time: %f\n", stop - start);

  check = check_matrix_array_equal(h_C_host, h_C, Coffsets, Coffsets, Ms_array, Ns_array, nBlocks);
  if( check == 0 ){
    printf("check host dgemm failed: %d\n",check);
  }

  //////////////////////////////////////////////////////////////////////
  // Run tests on device
  //////////////////////////////////////////////////////////////////////

  rocblas_handle handle;
  rocblas_create_handle(&handle);

  start = omp_get_wtime();
  rocblas_array_dgemm(handle, h_A, h_B, h_C_device, Ms_array, Ns_array, Ks_array, Aoffsets, Boffsets, Coffsets, nBlocks);
  stop = omp_get_wtime();
  printf("rocblas wrapper time: %f\n", stop - start);

  check = check_matrix_array_equal(h_C_device, h_C, Coffsets, Coffsets, Ms_array, Ns_array, nBlocks);
  if( check == 0 ){
    printf("check device dgemm failed: %d\n",check);
  }

  start = omp_get_wtime();
  magma_array_dgemm(h_A, h_B, h_C_magma, Ms_array, Ns_array, Ks_array, Aoffsets, Boffsets, Coffsets, nBlocks);
  stop = omp_get_wtime();
  printf("magma blas wrapper time: %f\n", stop - start);

  check = check_matrix_array_equal(h_C_magma, h_C, Coffsets, Coffsets, Ms_array, Ns_array, nBlocks);
  if( check == 0 ){
    printf("magma dgemm failed: %d\n",check);
  }

  // Clean up
  rocblas_destroy_handle(handle);
  free(h_A);
  free(h_B);
  free(h_C);
  free(h_C_host);
  free(h_C_device);
  free(h_C_magma);
  magma_finalize();

  return 0;
}

ext: c


text:
/***
 * Copyright 2019 The Katla Authors
 *
 * 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 KATLA_CORE_H
#define KATLA_CORE_H

#include "outcome/outcome.hpp"
namespace outcome = OUTCOME_V2_NAMESPACE;

#include "fmt/format.h"

#include <gsl/span>

#include <cstdio>
#include <cstdlib>
#include <string_view>

#ifdef _MSC_VER
    #ifdef KATLA_CORE_INDLL
    # define KATLA_CORE_DECLSPEC __declspec(dllexport)
    #else
    # define KATLA_CORE_DECLSPEC __declspec(dllimport)
    #endif
#else
    #define KATLA_CORE_DECLSPEC
#endif

namespace katla {

    template<class T, std::size_t Extent>
    using span = gsl::span<T, Extent>;

    /****
     * Declare format and print here for convenience. Current implementation uses the Fmt lib,
     * the idea is to use std c++ format later. Which is the reason for the indirection.
     */
    template <typename S, typename... Args>
    inline std::string format(const S& format_str, Args&&... args) {
        return fmt::format(format_str, args...);
    }

    template <typename S, typename... Args>
    inline void print(std::FILE* f, const S& format_str, Args&&... args) {
        fmt::print(f, format_str, args...);
    }

    template <typename S, typename... Args>
    inline void printInfo(const S& message, Args&&... args) {
        print(stdout, fmt::format("{}\n", message), args...);
        fflush(stdout);
    }

    template <typename S, typename... Args>
    inline void printError(const S& message, Args&&... args) {
        print(stderr, fmt::format("{}\n", message), args...);
        fflush(stderr);
    }

    template <typename S, typename... Args>
    inline void fatal(const S& message, Args&&... args) {
        print(stderr, fmt::format("{}\n", message), args...);
        fflush(stderr);
        std::abort();
    }
}

#endif // KATLA_CORE_H

max_stars_repo_name: plok/katla

max_issues_repo_name: plok/katla

max_forks_repo_name: plok/katla


text:
static char help[] = "Newton's method for a three-variable system.\n"
"The system is from fitting logistic population model\n"
"P(t) = N / (1 + A exp(-Nbt)) to 3 points (t,P) from census data in\n"
"D. Zill, 11th ed., exercise 4 in section 3.2.  The points are\n"
"(1790,3.929), (1850,23.192), (1910,91.972); populations in millions.\n";

/*
example:
$ make census
$ ./census -snes_fd -snes_monitor
  0 SNES Function norm 2.233793904584e+01
...
 17 SNES Function norm 6.521724317021e-07
 18 SNES Function norm 4.358557709045e-09
Vec Object: 1 MPI processes
  type: seq
197.274
49.2096
0.000158863
*/

#include <petsc.h>

extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    SNES  snes;          // nonlinear solver context
    Vec   x, r;          // solution, residual vectors
    double *ax;

    PetscInitialize(&argc,&argv,NULL,help);
    ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,3); CHKERRQ(ierr);
    ierr = VecSetFromOptions(x); CHKERRQ(ierr);
    ierr = VecSetUp(x); CHKERRQ(ierr);
    ierr = VecGetArray(x,&ax);CHKERRQ(ierr);
    // ax = [N, A, b]; see equations in FormFunction() below
    // these initial values are based on *guessing* N=300
    // as the limiting population, then computing A, then b
    ax[0] = 300.0;
    ax[1] = 75.0;
    ax[2] = 1.0e-4;
    ierr = VecRestoreArray(x,&ax);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(x);CHKERRQ(ierr);
    ierr = VecDuplicate(x,&r); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,r,FormFunction,NULL); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,x); CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    VecDestroy(&x);  VecDestroy(&r);  SNESDestroy(&snes);
    return PetscFinalize();
}

PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void *ctx) {
    PetscErrorCode ierr;
    const double *ax;
    double       *aF;

    ierr = VecGetArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF);CHKERRQ(ierr);
    // ax = [N, A, b]
    //  3.929 (1 + A)               - N = 0
    // 23.192 (1 + A exp(-60*N*b))  - N = 0
    // 91.972 (1 + A exp(-120*N*b)) - N = 0
    aF[0] =  3.929 * (1.0 + ax[1]) - ax[0];
    aF[1] = 23.192 * (1.0 + ax[1] * PetscExpReal(-60.0*ax[0]*ax[2])) - ax[0];
    aF[2] = 91.972 * (1.0 + ax[1] * PetscExpReal(-120.0*ax[0]*ax[2])) - ax[0];
    ierr = VecRestoreArray(F,&aF);CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(x,&ax);CHKERRQ(ierr);
    return 0;
}


text:
/* multimin/simplex.c
 *
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 2002 Tuomo Keskitalo, Ivo Alxneit
 *
 * 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.
 */

/*
   - Originally written by Tuomo Keskitalo <tuomo.keskitalo@iki.fi>
   - Corrections to nmsimplex_iterate and other functions
     by Ivo Alxneit <ivo.alxneit@psi.ch>
   - Additional help by Brian Gough <bjg@network-theory.co.uk>
*/

/* The Simplex method of Nelder and Mead,
   also known as the polytope search alogorithm. Ref:
   Nelder, J.A., Mead, R., Computer Journal 7 (1965) pp. 308-313.

   This implementation uses n+1 corner points in the simplex.
*/

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multimin.h>

typedef struct
{
  gsl_matrix *x1;               /* simplex corner points */
  gsl_vector *y1;               /* function value at corner points */
  gsl_vector *ws1;              /* workspace 1 for algorithm */
  gsl_vector *ws2;              /* workspace 2 for algorithm */
}
nmsimplex_state_t;

static double
nmsimplex_move_corner (const double coeff, const nmsimplex_state_t * state,
                       size_t corner, gsl_vector * xc,
                       const gsl_multimin_function * f)
{
  /* moves a simplex corner scaled by coeff (negative value represents
     mirroring by the middle point of the "other" corner points)
     and gives new corner in xc and function value at xc as a
     return value
   */

  gsl_matrix *x1 = state->x1;

  size_t i, j;
  double newval, mp;

  for (j = 0; j < x1->size2; j++)
    {
      mp = 0.0;
      for (i = 0; i < x1->size1; i++)
        {
          if (i != corner)
            {
              mp += (gsl_matrix_get (x1, i, j));
            }
        }
      mp /= (double) (x1->size1 - 1);
      newval = mp - coeff * (mp - gsl_matrix_get (x1, corner, j));
      gsl_vector_set (xc, j, newval);
    }

  newval = GSL_MULTIMIN_FN_EVAL (f, xc);

  return newval;
}

static int
nmsimplex_contract_by_best (nmsimplex_state_t * state, size_t best,
                            gsl_vector * xc, gsl_multimin_function * f)
{

  /* Function contracts the simplex in respect to
     best valued corner. That is, all corners besides the
     best corner are moved. */

  /* the xc vector is simply work space here */

  gsl_matrix *x1 = state->x1;
  gsl_vector *y1 = state->y1;

  size_t i, j;
  double newval;

  int status = GSL_SUCCESS;

  for (i = 0; i < x1->size1; i++)
    {
      if (i != best)
        {
          for (j = 0; j < x1->size2; j++)
            {
              newval = 0.5 * (gsl_matrix_get (x1, i, j)
                              + gsl_matrix_get (x1, best, j));
              gsl_matrix_set (x1, i, j, newval);
            }

          /* evaluate function in the new point */

          gsl_matrix_get_row (xc, x1, i);
          newval = GSL_MULTIMIN_FN_EVAL (f, xc);
          gsl_vector_set (y1, i, newval);

	  /* notify caller that we found at least one bad function value.
	     we finish the contraction (and do not abort) to allow the user
	     to handle the situation */

          if(!gsl_finite(newval))
	    {
	      status = GSL_EBADFUNC;
	    }
        }
    }

  return status;
}

static int
nmsimplex_calc_center (const nmsimplex_state_t * state, gsl_vector * mp)
{
  /* calculates the center of the simplex to mp */

  gsl_matrix *x1 = state->x1;

  size_t i, j;
  double val;

  for (j = 0; j < x1->size2; j++)
    {
      val = 0.0;
      for (i = 0; i < x1->size1; i++)
        {
          val += gsl_matrix_get (x1, i, j);
        }
      val /= x1->size1;
      gsl_vector_set (mp, j, val);
    }

  return GSL_SUCCESS;
}

static double
nmsimplex_size (nmsimplex_state_t * state)
{
  /* calculates simplex size as average sum of length of vectors
     from simplex center to corner points:

     ( sum ( || y - y_middlepoint || ) ) / n
   */

  gsl_vector *s = state->ws1;
  gsl_vector *mp = state->ws2;

  gsl_matrix *x1 = state->x1;
  size_t i;

  double ss = 0.0;

  /* Calculate middle point */
  nmsimplex_calc_center (state, mp);

  for (i = 0; i < x1->size1; i++)
    {
      gsl_matrix_get_row (s, x1, i);
      gsl_blas_daxpy (-1.0, mp, s);
      ss += gsl_blas_dnrm2 (s);
    }

  return ss / (double) (x1->size1);
}

static int
nmsimplex_alloc (void *vstate, size_t n)
{
  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  if (n == 0)
    {
      GSL_ERROR("invalid number of parameters specified", GSL_EINVAL);
    }

  state->x1 = gsl_matrix_alloc (n + 1, n);

  if (state->x1 == NULL)
    {
      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
    }

  state->y1 = gsl_vector_alloc (n + 1);

  if (state->y1 == NULL)
    {
      gsl_matrix_free(state->x1);
      GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
    }

  state->ws1 = gsl_vector_alloc (n);

  if (state->ws1 == NULL)
    {
      gsl_matrix_free(state->x1);
      gsl_vector_free(state->y1);
      GSL_ERROR ("failed to allocate space for ws1", GSL_ENOMEM);
    }

  state->ws2 = gsl_vector_alloc (n);

  if (state->ws2 == NULL)
    {
      gsl_matrix_free(state->x1);
      gsl_vector_free(state->y1);
      gsl_vector_free(state->ws1);
      GSL_ERROR ("failed to allocate space for ws2", GSL_ENOMEM);
    }

  return GSL_SUCCESS;
}

static int
nmsimplex_set (void *vstate, gsl_multimin_function * f,
               const gsl_vector * x,
               double *size, const gsl_vector * step_size)
{
  int status;
  size_t i;
  double val;

  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  gsl_vector *xtemp = state->ws1;

  if (xtemp->size != x->size)
    {
      GSL_ERROR("incompatible size of x", GSL_EINVAL);
    }

  if (xtemp->size != step_size->size)
    {
      GSL_ERROR("incompatible size of step_size", GSL_EINVAL);
    }

  /* first point is the original x0 */

  val = GSL_MULTIMIN_FN_EVAL (f, x);

  if (!gsl_finite(val))
    {
      GSL_ERROR("non-finite function value encountered", GSL_EBADFUNC);
    }

  gsl_matrix_set_row (state->x1, 0, x);
  gsl_vector_set (state->y1, 0, val);

  /* following points are initialized to x0 + step_size */

  for (i = 0; i < x->size; i++)
    {
      status = gsl_vector_memcpy (xtemp, x);

      if (status != 0)
        {
          GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
        }

      val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i);
      gsl_vector_set (xtemp, i, val);
      val = GSL_MULTIMIN_FN_EVAL (f, xtemp);

      if (!gsl_finite(val))
        {
          GSL_ERROR("non-finite function value encountered", GSL_EBADFUNC);
        }

      gsl_matrix_set_row (state->x1, i + 1, xtemp);
      gsl_vector_set (state->y1, i + 1, val);
    }

  /* Initialize simplex size */

  *size = nmsimplex_size (state);

  return GSL_SUCCESS;
}

static void
nmsimplex_free (void *vstate)
{
  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  gsl_matrix_free (state->x1);
  gsl_vector_free (state->y1);
  gsl_vector_free (state->ws1);
  gsl_vector_free (state->ws2);
}

static int
nmsimplex_iterate (void *vstate, gsl_multimin_function * f,
                   gsl_vector * x, double *size, double *fval)
{

  /* Simplex iteration tries to minimize function f value */
  /* Includes corrections from Ivo Alxneit <ivo.alxneit@psi.ch> */

  nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;

  /* xc and xc2 vectors store tried corner point coordinates */

  gsl_vector *xc = state->ws1;
  gsl_vector *xc2 = state->ws2;
  gsl_vector *y1 = state->y1;
  gsl_matrix *x1 = state->x1;

  size_t n = y1->size;
  size_t i;
  size_t hi, s_hi, lo;
  double dhi, ds_hi, dlo;
  int status;
  double val, val2;

  if (xc->size != x->size)
    {
      GSL_ERROR("incompatible size of x", GSL_EINVAL);
    }

  /* get index of highest, second highest and lowest point */

  dhi = dlo = gsl_vector_get (y1, 0);
  hi = 0; lo = 0;

  ds_hi = gsl_vector_get(y1, 1);
  s_hi = 1;

  for (i = 1; i < n; i++)
    {
      val = (gsl_vector_get (y1, i));
      if (val < dlo)
        {
          dlo = val;
          lo = i;
        }
      else if (val > dhi)
        {
          ds_hi = dhi;
          s_hi = hi;
          dhi = val;
          hi = i;
        }
      else if (val > ds_hi)
        {
          ds_hi = val;
          s_hi = i;
        }
    }

  /* reflect the highest value */

  val = nmsimplex_move_corner (-1.0, state, hi, xc, f);

  if (gsl_finite(val) && val < gsl_vector_get (y1, lo))
    {

      /* reflected point becomes lowest point, try expansion */

      val2 = nmsimplex_move_corner (-2.0, state, hi, xc2, f);

      if (gsl_finite(val2) && val2 < gsl_vector_get (y1, lo))
        {
          gsl_matrix_set_row (x1, hi, xc2);
          gsl_vector_set (y1, hi, val2);
        }
      else
        {
          gsl_matrix_set_row (x1, hi, xc);
          gsl_vector_set (y1, hi, val);
        }
    }

  /* reflection does not improve things enough
     or
     we got a non-finite (illegal) function value */

  else if (!gsl_finite(val) || val > gsl_vector_get (y1, s_hi))
    {
      if (gsl_finite(val) && val <= gsl_vector_get (y1, hi))
        {

          /* if trial point is better than highest point, replace
             highest point */

          gsl_matrix_set_row (x1, hi, xc);
          gsl_vector_set (y1, hi, val);
        }

      /* try one dimensional contraction */

      val2 = nmsimplex_move_corner (0.5, state, hi, xc2, f);

      if (gsl_finite(val2) && val2 <= gsl_vector_get (y1, hi))
        {
          gsl_matrix_set_row (state->x1, hi, xc2);
          gsl_vector_set (y1, hi, val2);
        }

      else
        {

          /* contract the whole simplex in respect to the best point */

          status = nmsimplex_contract_by_best (state, lo, xc, f);
          if (status != GSL_SUCCESS)
            {
              GSL_ERROR ("nmsimplex_contract_by_best failed", GSL_EFAILED);
            }
        }
    }
  else
    {

      /* trial point is better than second highest point.
         Replace highest point by it */

      gsl_matrix_set_row (x1, hi, xc);
      gsl_vector_set (y1, hi, val);
    }

  /* return lowest point of simplex as x */

  lo = gsl_vector_min_index (y1);
  gsl_matrix_get_row (x, x1, lo);
  *fval = gsl_vector_get (y1, lo);

  /* Update simplex size */

  *size = nmsimplex_size (state);

  return GSL_SUCCESS;
}

static const gsl_multimin_fminimizer_type nmsimplex_type =
{ "nmsimplex",  /* name */
  sizeof (nmsimplex_state_t),
  &nmsimplex_alloc,
  &nmsimplex_set,
  &nmsimplex_iterate,
  &nmsimplex_free
};

const gsl_multimin_fminimizer_type
  * gsl_multimin_fminimizer_nmsimplex = &nmsimplex_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
//        Copyright The Authors 2018.
//    Distributed under the 3-Clause BSD License.
//    (See accompanying file LICENSE or copy at
//   https://opensource.org/licenses/BSD-3-Clause)

#pragma once

#include <cstdint>  // for std::uint8_t
#include <gsl/gsl>  // for gsl::span<std::uint8_t>
#include <vector>   // for std::vector<std::uint8_t>

namespace blocxxi {
namespace p2p {
namespace kademlia {

/// Represents an output buffer to which headers and message bodies can be
/// serialized. Whatever type is used, it needs to be able to expand as needed.
using Buffer = std::vector<std::uint8_t>;

/// Represents a read-only view over the input buffer used to deserialize
/// headers and message bodies.
using BufferReader = gsl::span<std::uint8_t>;

}  // namespace kademlia
}  // namespace p2p
}  // namespace blocxxi


text:
#include <stdio.h>
#include <gsl/gsl_const_mksa.h>

int
main (void)
{
  double c  = GSL_CONST_MKS_SPEED_OF_LIGHT;
  double au = GSL_CONST_MKS_ASTRONOMICAL_UNIT;
  double min = GSL_CONST_MKS_MINUTE;

  double r_earth = 1.00 * au;
  double r_mars  = 1.52 * au;

  double tmin, tmax;

  tmin = (r_mars - r_earth)/c;
  tmax = (r_mars + r_earth)/c;

  printf("light travel time from Earth to Mars:\n");
  printf("minimum = %.1f minutes\n", tmin/min);
  printf("maximum = %.1f minutes\n", tmax/min);
}

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:
/* matrix/gsl_matrix_ushort.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 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_MATRIX_USHORT_H__
#define __GSL_MATRIX_USHORT_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_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_ushort.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;
  unsigned short * data;
  gsl_block_ushort * block;
  int owner;
} gsl_matrix_ushort;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_view;

typedef _gsl_matrix_ushort_view gsl_matrix_ushort_view;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_const_view;

typedef const _gsl_matrix_ushort_const_view gsl_matrix_ushort_const_view;

/* Allocation */

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_block (gsl_block_ushort * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_matrix (gsl_matrix_ushort * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_ushort *
gsl_vector_ushort_alloc_row_from_matrix (gsl_matrix_ushort * m,
                                        const size_t i);

GSL_FUN gsl_vector_ushort *
gsl_vector_ushort_alloc_col_from_matrix (gsl_matrix_ushort * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_ushort_free (gsl_matrix_ushort * m);

/* Views */

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_submatrix (gsl_matrix_ushort * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_row (gsl_matrix_ushort * m, const size_t i);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_column (gsl_matrix_ushort * m, const size_t j);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_diagonal (gsl_matrix_ushort * m);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_subdiagonal (gsl_matrix_ushort * m, const size_t k);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_superdiagonal (gsl_matrix_ushort * m, const size_t k);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_subrow (gsl_matrix_ushort * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_subcolumn (gsl_matrix_ushort * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_array (unsigned short * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_array_with_tda (unsigned short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector (gsl_vector_ushort * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector_with_tda (gsl_vector_ushort * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_submatrix (const gsl_matrix_ushort * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_row (const gsl_matrix_ushort * m,
                            const size_t i);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_column (const gsl_matrix_ushort * m,
                               const size_t j);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_diagonal (const gsl_matrix_ushort * m);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subdiagonal (const gsl_matrix_ushort * m,
                                    const size_t k);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_superdiagonal (const gsl_matrix_ushort * m,
                                      const size_t k);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subrow (const gsl_matrix_ushort * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subcolumn (const gsl_matrix_ushort * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array (const unsigned short * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array_with_tda (const unsigned short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector (const gsl_vector_ushort * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector_with_tda (const gsl_vector_ushort * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_ushort_set_zero (gsl_matrix_ushort * m);
GSL_FUN void gsl_matrix_ushort_set_identity (gsl_matrix_ushort * m);
GSL_FUN void gsl_matrix_ushort_set_all (gsl_matrix_ushort * m, unsigned short x);

GSL_FUN int gsl_matrix_ushort_fread (FILE * stream, gsl_matrix_ushort * m) ;
GSL_FUN int gsl_matrix_ushort_fwrite (FILE * stream, const gsl_matrix_ushort * m) ;
GSL_FUN int gsl_matrix_ushort_fscanf (FILE * stream, gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_fprintf (FILE * stream, const gsl_matrix_ushort * m, const char * format);

GSL_FUN int gsl_matrix_ushort_memcpy(gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);
GSL_FUN int gsl_matrix_ushort_swap(gsl_matrix_ushort * m1, gsl_matrix_ushort * m2);

GSL_FUN int gsl_matrix_ushort_swap_rows(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ushort_swap_columns(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ushort_swap_rowcol(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ushort_transpose (gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_transpose_memcpy (gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);

GSL_FUN unsigned short gsl_matrix_ushort_max (const gsl_matrix_ushort * m);
GSL_FUN unsigned short gsl_matrix_ushort_min (const gsl_matrix_ushort * m);
GSL_FUN void gsl_matrix_ushort_minmax (const gsl_matrix_ushort * m, unsigned short * min_out, unsigned short * max_out);

GSL_FUN void gsl_matrix_ushort_max_index (const gsl_matrix_ushort * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_ushort_min_index (const gsl_matrix_ushort * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_ushort_minmax_index (const gsl_matrix_ushort * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_ushort_equal (const gsl_matrix_ushort * a, const gsl_matrix_ushort * b);

GSL_FUN int gsl_matrix_ushort_isnull (const gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_ispos (const gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_isneg (const gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_isnonneg (const gsl_matrix_ushort * m);

GSL_FUN int gsl_matrix_ushort_add (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_sub (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_mul_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_div_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_scale (gsl_matrix_ushort * a, const double x);
GSL_FUN int gsl_matrix_ushort_add_constant (gsl_matrix_ushort * a, const double x);
GSL_FUN int gsl_matrix_ushort_add_diagonal (gsl_matrix_ushort * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_ushort_get_row(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t i);
GSL_FUN int gsl_matrix_ushort_get_col(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t j);
GSL_FUN int gsl_matrix_ushort_set_row(gsl_matrix_ushort * m, const size_t i, const gsl_vector_ushort * v);
GSL_FUN int gsl_matrix_ushort_set_col(gsl_matrix_ushort * m, const size_t j, const gsl_vector_ushort * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL unsigned short   gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x);
GSL_FUN INLINE_DECL unsigned short * gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const unsigned short * gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
unsigned short
gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      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
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
unsigned short *
gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      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 (unsigned short *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const unsigned short *
gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      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 unsigned short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_USHORT_H__ */


text:
#pragma once

#include <GLES2/gl2.h>
#include <GLES2/gl2ext.h>
#include <GLES3/gl3.h>
#include <GLES3/gl3ext.h>
#include <GLES3/gl3platform.h>
#include <EGL/egl.h>
#include <EGL/eglext.h>
#include <gsl/gsl>

namespace android::OpenGLHelpers
{
    constexpr GLint GetTextureUnit(GLenum texture)
    {
        return texture - GL_TEXTURE0;
    }

    GLuint CreateShaderProgram(const char* vertShaderSource, const char* fragShaderSource);

    namespace GLTransactions
    {
        inline auto MakeCurrent(EGLDisplay display, EGLSurface drawSurface, EGLSurface readSurface, EGLContext context)
        {
            EGLDisplay previousDisplay{ eglGetDisplay(EGL_DEFAULT_DISPLAY) };
            EGLSurface previousDrawSurface{ eglGetCurrentSurface(EGL_DRAW) };
            EGLSurface previousReadSurface{ eglGetCurrentSurface(EGL_READ) };
            EGLContext previousContext{ eglGetCurrentContext() };
            eglMakeCurrent(display, drawSurface, readSurface, context);
            return gsl::finally([previousDisplay, previousDrawSurface, previousReadSurface, previousContext]() { eglMakeCurrent(previousDisplay, previousDrawSurface, previousReadSurface, previousContext); });
        }
    }
}


text:
/* linalg/svd.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 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 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_blas.h>

#include <gsl/gsl_linalg.h>

#include "givens.c"
#include "svdstep.c"

/* Factorise a general M x N matrix A into,
 *
 *   A = U D V^T
 *
 * where U is a column-orthogonal M x N matrix (U^T U = I),
 * D is a diagonal N x N matrix,
 * and V is an N x N orthogonal matrix (V^T V = V V^T = I)
 *
 * U is stored in the original matrix A, which has the same size
 *
 * V is stored as a separate matrix (not V^T). You must take the
 * transpose to form the product above.
 *
 * The diagonal matrix D is stored in the vector S,  D_ii = S_i
 */

int
gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S,
                      gsl_vector * work)
{
  size_t a, b, i, j, iter;

  const size_t M = A->size1;
  const size_t N = A->size2;
  const size_t K = GSL_MIN (M, N);

  if (M < N)
    {
      GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
    }
  else if (V->size1 != N)
    {
      GSL_ERROR ("square matrix V must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else if (V->size1 != V->size2)
    {
      GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
    }
  else if (S->size != N)
    {
      GSL_ERROR ("length of vector S must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else if (work->size != N)
    {
      GSL_ERROR ("length of workspace must match second dimension of matrix A",
                 GSL_EBADLEN);
    }

  /* Handle the case of N = 1 (SVD of a column vector) */

  if (N == 1)
    {
      gsl_vector_view column = gsl_matrix_column (A, 0);
      double norm = gsl_blas_dnrm2 (&column.vector);

      gsl_vector_set (S, 0, norm);
      gsl_matrix_set (V, 0, 0, 1.0);

      if (norm != 0.0)
        {
          gsl_blas_dscal (1.0/norm, &column.vector);
        }

      return GSL_SUCCESS;
    }

  {
    gsl_vector_view f = gsl_vector_subvector (work, 0, K - 1);

    /* bidiagonalize matrix A, unpack A into U S V */

    gsl_linalg_bidiag_decomp (A, S, &f.vector);
    gsl_linalg_bidiag_unpack2 (A, S, &f.vector, V);

    /* apply reduction steps to B=(S,Sd) */

    chop_small_elements (S, &f.vector);

    /* Progressively reduce the matrix until it is diagonal */

    b = N - 1;
    iter = 0;

    while (b > 0)
      {
        double fbm1 = gsl_vector_get (&f.vector, b - 1);

        if (fbm1 == 0.0 || gsl_isnan (fbm1))
          {
            b--;
            continue;
          }

        /* Find the largest unreduced block (a,b) starting from b
           and working backwards */

        a = b - 1;

        while (a > 0)
          {
            double fam1 = gsl_vector_get (&f.vector, a - 1);

            if (fam1 == 0.0 || gsl_isnan (fam1))
              {
                break;
              }

            a--;
          }

        iter++;

        if (iter > 100 * N)
          {
            GSL_ERROR("SVD decomposition failed to converge", GSL_EMAXITER);
          }

        {
          const size_t n_block = b - a + 1;
          gsl_vector_view S_block = gsl_vector_subvector (S, a, n_block);
          gsl_vector_view f_block = gsl_vector_subvector (&f.vector, a, n_block - 1);

          gsl_matrix_view U_block =
            gsl_matrix_submatrix (A, 0, a, A->size1, n_block);
          gsl_matrix_view V_block =
            gsl_matrix_submatrix (V, 0, a, V->size1, n_block);

          qrstep (&S_block.vector, &f_block.vector, &U_block.matrix, &V_block.matrix);

          /* remove any small off-diagonal elements */

          chop_small_elements (&S_block.vector, &f_block.vector);
        }
      }
  }
  /* Make singular values positive by reflections if necessary */

  for (j = 0; j < K; j++)
    {
      double Sj = gsl_vector_get (S, j);

      if (Sj < 0.0)
        {
          for (i = 0; i < N; i++)
            {
              double Vij = gsl_matrix_get (V, i, j);
              gsl_matrix_set (V, i, j, -Vij);
            }

          gsl_vector_set (S, j, -Sj);
        }
    }

  /* Sort singular values into decreasing order */

  for (i = 0; i < K; i++)
    {
      double S_max = gsl_vector_get (S, i);
      size_t i_max = i;

      for (j = i + 1; j < K; j++)
        {
          double Sj = gsl_vector_get (S, j);

          if (Sj > S_max)
            {
              S_max = Sj;
              i_max = j;
            }
        }

      if (i_max != i)
        {
          /* swap eigenvalues */
          gsl_vector_swap_elements (S, i, i_max);

          /* swap eigenvectors */
          gsl_matrix_swap_columns (A, i, i_max);
          gsl_matrix_swap_columns (V, i, i_max);
        }
    }

  return GSL_SUCCESS;
}

/* Modified algorithm which is better for M>>N */

int
gsl_linalg_SV_decomp_mod (gsl_matrix * A,
                          gsl_matrix * X,
                          gsl_matrix * V, gsl_vector * S, gsl_vector * work)
{
  size_t i, j;

  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M < N)
    {
      GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
    }
  else if (V->size1 != N)
    {
      GSL_ERROR ("square matrix V must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else if (V->size1 != V->size2)
    {
      GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
    }
  else if (X->size1 != N)
    {
      GSL_ERROR ("square matrix X must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else if (X->size1 != X->size2)
    {
      GSL_ERROR ("matrix X must be square", GSL_ENOTSQR);
    }
  else if (S->size != N)
    {
      GSL_ERROR ("length of vector S must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else if (work->size != N)
    {
      GSL_ERROR ("length of workspace must match second dimension of matrix A",
                 GSL_EBADLEN);
    }

  if (N == 1)
    {
      gsl_vector_view column = gsl_matrix_column (A, 0);
      double norm = gsl_blas_dnrm2 (&column.vector);

      gsl_vector_set (S, 0, norm);
      gsl_matrix_set (V, 0, 0, 1.0);

      if (norm != 0.0)
        {
          gsl_blas_dscal (1.0/norm, &column.vector);
        }

      return GSL_SUCCESS;
    }

  /* Convert A into an upper triangular matrix R */

  for (i = 0; i < N; i++)
    {
      gsl_vector_view c = gsl_matrix_column (A, i);
      gsl_vector_view v = gsl_vector_subvector (&c.vector, i, M - i);
      double tau_i = gsl_linalg_householder_transform (&v.vector);

      /* Apply the transformation to the remaining columns */

      if (i + 1 < N)
        {
          gsl_matrix_view m =
            gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
          gsl_linalg_householder_hm (tau_i, &v.vector, &m.matrix);
        }

      gsl_vector_set (S, i, tau_i);
    }

  /* Copy the upper triangular part of A into X */

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < i; j++)
        {
          gsl_matrix_set (X, i, j, 0.0);
        }

      {
        double Aii = gsl_matrix_get (A, i, i);
        gsl_matrix_set (X, i, i, Aii);
      }

      for (j = i + 1; j < N; j++)
        {
          double Aij = gsl_matrix_get (A, i, j);
          gsl_matrix_set (X, i, j, Aij);
        }
    }

  /* Convert A into an orthogonal matrix L */

  for (j = N; j-- > 0;)
    {
      /* Householder column transformation to accumulate L */
      double tj = gsl_vector_get (S, j);
      gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M - j, N - j);
      gsl_linalg_householder_hm1 (tj, &m.matrix);
    }

  /* unpack R into X V S */

  gsl_linalg_SV_decomp (X, V, S, work);

  /* Multiply L by X, to obtain U = L X, stored in U */

  {
    gsl_vector_view sum = gsl_vector_subvector (work, 0, N);

    for (i = 0; i < M; i++)
      {
        gsl_vector_view L_i = gsl_matrix_row (A, i);
        gsl_vector_set_zero (&sum.vector);

        for (j = 0; j < N; j++)
          {
            double Lij = gsl_vector_get (&L_i.vector, j);
            gsl_vector_view X_j = gsl_matrix_row (X, j);
            gsl_blas_daxpy (Lij, &X_j.vector, &sum.vector);
          }

        gsl_vector_memcpy (&L_i.vector, &sum.vector);
      }
  }

  return GSL_SUCCESS;
}

/*  Solves the system A x = b using the SVD factorization
 *
 *  A = U S V^T
 *
 *  to obtain x. For M x N systems it finds the solution in the least
 *  squares sense.
 */

int
gsl_linalg_SV_solve (const gsl_matrix * U,
                     const gsl_matrix * V,
                     const gsl_vector * S,
                     const gsl_vector * b, gsl_vector * x)
{
  if (U->size1 != b->size)
    {
      GSL_ERROR ("first dimension of matrix U must size of vector b",
                 GSL_EBADLEN);
    }
  else if (U->size2 != S->size)
    {
      GSL_ERROR ("length of vector S must match second dimension of matrix U",
                 GSL_EBADLEN);
    }
  else if (V->size1 != V->size2)
    {
      GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
    }
  else if (S->size != V->size1)
    {
      GSL_ERROR ("length of vector S must match size of matrix V",
                 GSL_EBADLEN);
    }
  else if (V->size2 != x->size)
    {
      GSL_ERROR ("size of matrix V must match size of vector x", GSL_EBADLEN);
    }
  else
    {
      const size_t N = U->size2;
      size_t i;

      gsl_vector *w = gsl_vector_calloc (N);

      gsl_blas_dgemv (CblasTrans, 1.0, U, b, 0.0, w);

      for (i = 0; i < N; i++)
        {
          double wi = gsl_vector_get (w, i);
          double alpha = gsl_vector_get (S, i);
          if (alpha != 0)
            alpha = 1.0 / alpha;
          gsl_vector_set (w, i, alpha * wi);
        }

      gsl_blas_dgemv (CblasNoTrans, 1.0, V, w, 0.0, x);

      gsl_vector_free (w);

      return GSL_SUCCESS;
    }
}

/* This is a the jacobi version */
/* Author:  G. Jungman */

/*
 * Algorithm due to J.C. Nash, Compact Numerical Methods for
 * Computers (New York: Wiley and Sons, 1979), chapter 3.
 * See also Algorithm 4.1 in
 * James Demmel, Kresimir Veselic, "Jacobi's Method is more
 * accurate than QR", Lapack Working Note 15 (LAWN15), October 1989.
 * Available from netlib.
 *
 * Based on code by Arthur Kosowsky, Rutgers University
 *                  kosowsky@physics.rutgers.edu
 *
 * Another relevant paper is, P.P.M. De Rijk, "A One-Sided Jacobi
 * Algorithm for computing the singular value decomposition on a
 * vector computer", SIAM Journal of Scientific and Statistical
 * Computing, Vol 10, No 2, pp 359-371, March 1989.
 *
 */

int
gsl_linalg_SV_decomp_jacobi (gsl_matrix * A, gsl_matrix * Q, gsl_vector * S)
{
  if (A->size1 < A->size2)
    {
      /* FIXME: only implemented  M>=N case so far */

      GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
    }
  else if (Q->size1 != A->size2)
    {
      GSL_ERROR ("square matrix Q must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else if (Q->size1 != Q->size2)
    {
      GSL_ERROR ("matrix Q must be square", GSL_ENOTSQR);
    }
  else if (S->size != A->size2)
    {
      GSL_ERROR ("length of vector S must match second dimension of matrix A",
                 GSL_EBADLEN);
    }
  else
    {
      const size_t M = A->size1;
      const size_t N = A->size2;
      size_t i, j, k;

      /* Initialize the rotation counter and the sweep counter. */
      int count = 1;
      int sweep = 0;
      int sweepmax = 5*N;

      double tolerance = 10 * M * GSL_DBL_EPSILON;

      /* Always do at least 12 sweeps. */
      sweepmax = GSL_MAX (sweepmax, 12);

      /* Set Q to the identity matrix. */
      gsl_matrix_set_identity (Q);

      /* Store the column error estimates in S, for use during the
         orthogonalization */

      for (j = 0; j < N; j++)
        {
          gsl_vector_view cj = gsl_matrix_column (A, j);
          double sj = gsl_blas_dnrm2 (&cj.vector);
          gsl_vector_set(S, j, GSL_DBL_EPSILON * sj);
        }

      /* Orthogonalize A by plane rotations. */

      while (count > 0 && sweep <= sweepmax)
        {
          /* Initialize rotation counter. */
          count = N * (N - 1) / 2;

          for (j = 0; j < N - 1; j++)
            {
              for (k = j + 1; k < N; k++)
                {
                  double a = 0.0;
                  double b = 0.0;
                  double p = 0.0;
                  double q = 0.0;
                  double cosine, sine;
                  double v;
                  double abserr_a, abserr_b;
                  int sorted, orthog, noisya, noisyb;

                  gsl_vector_view cj = gsl_matrix_column (A, j);
                  gsl_vector_view ck = gsl_matrix_column (A, k);

                  gsl_blas_ddot (&cj.vector, &ck.vector, &p);
                  p *= 2.0 ;  /* equation 9a:  p = 2 x.y */

                  a = gsl_blas_dnrm2 (&cj.vector);
                  b = gsl_blas_dnrm2 (&ck.vector);

                  q = a * a - b * b;
                  v = hypot(p, q);

                  /* test for columns j,k orthogonal, or dominant errors */

                  abserr_a = gsl_vector_get(S,j);
                  abserr_b = gsl_vector_get(S,k);

                  sorted = (GSL_COERCE_DBL(a) >= GSL_COERCE_DBL(b));
                  orthog = (fabs (p) <= tolerance * GSL_COERCE_DBL(a * b));
                  noisya = (a < abserr_a);
                  noisyb = (b < abserr_b);

                  if (sorted && (orthog || noisya || noisyb))
                    {
                      count--;
                      continue;
                    }

                  /* calculate rotation angles */
                  if (v == 0 || !sorted)
                    {
                      cosine = 0.0;
                      sine = 1.0;
                    }
                  else
                    {
                      cosine = sqrt((v + q) / (2.0 * v));
                      sine = p / (2.0 * v * cosine);
                    }

                  /* apply rotation to A */
                  for (i = 0; i < M; i++)
                    {
                      const double Aik = gsl_matrix_get (A, i, k);
                      const double Aij = gsl_matrix_get (A, i, j);
                      gsl_matrix_set (A, i, j, Aij * cosine + Aik * sine);
                      gsl_matrix_set (A, i, k, -Aij * sine + Aik * cosine);
                    }

                  gsl_vector_set(S, j, fabs(cosine) * abserr_a + fabs(sine) * abserr_b);
                  gsl_vector_set(S, k, fabs(sine) * abserr_a + fabs(cosine) * abserr_b);

                  /* apply rotation to Q */
                  for (i = 0; i < N; i++)
                    {
                      const double Qij = gsl_matrix_get (Q, i, j);
                      const double Qik = gsl_matrix_get (Q, i, k);
                      gsl_matrix_set (Q, i, j, Qij * cosine + Qik * sine);
                      gsl_matrix_set (Q, i, k, -Qij * sine + Qik * cosine);
                    }
                }
            }

          /* Sweep completed. */
          sweep++;
        }

      /*
       * Orthogonalization complete. Compute singular values.
       */

      {
        double prev_norm = -1.0;

        for (j = 0; j < N; j++)
          {
            gsl_vector_view column = gsl_matrix_column (A, j);
            double norm = gsl_blas_dnrm2 (&column.vector);

            /* Determine if singular value is zero, according to the
               criteria used in the main loop above (i.e. comparison
               with norm of previous column). */

            if (norm == 0.0 || prev_norm == 0.0
                || (j > 0 && norm <= tolerance * prev_norm))
              {
                gsl_vector_set (S, j, 0.0);     /* singular */
                gsl_vector_set_zero (&column.vector);   /* annihilate column */

                prev_norm = 0.0;
              }
            else
              {
                gsl_vector_set (S, j, norm);    /* non-singular */
                gsl_vector_scale (&column.vector, 1.0 / norm);  /* normalize column */

                prev_norm = norm;
              }
          }
      }

      if (count > 0)
        {
          /* reached sweep limit */
          GSL_ERROR ("Jacobi iterations did not reach desired tolerance",
                     GSL_ETOL);
        }

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: mattjr/structured


text:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <math.h>
#ifndef CRANESUB
char *fn1="crane1io.dat";
char *fn2="crane2io.dat";
FILE *fp;
#endif//#ifndef CRANESUB
typedef struct{
  double g;//gravity
  double M;//mass of the car
  double m;//mass of the load
  double r;//length of the rope
  double C;//damping coefficient

  double F;//input force
  double T;//tension to the rope

  double x,x0;//x position of the load
  double y,y0;//y position of the load
  double X;//X position of the car
  double a;//angle (theta)

  double dx,dx0;
  double dy,dy0;
  double dX,dX0;
  double dr;
  double da,da0;

  double ddx;
  double ddy;
  double ddX;
  double ddr;
  double dda;
  double dXmax;
  double ddXmax;
  double Fmax;//Force

  double h;
  int nh;
  double _h;
  int _dim;
  double *_y;
  double *_y_err;
  double *_dydt_in;
  double *_dydt_out;
  const gsl_odeiv_step_type *_T;
  gsl_odeiv_step *_s;
  gsl_odeiv_system _sys;
  double _t;
  //  int *initialize;
  int (*cranefunc)();
  double (*plant)();
  //int (*p) (int num);    /* 関数ポインタpを宣言 */
} CRANE;

#define dim_crane 4
//double ddXn;//for check
int crane1func (double t, const double y[], double f[], void *params)
{
  CRANE *c = (CRANE *)params;
  //y[0]=crane->p, y[1]=crane->dp, y[2]=crane->X, y[3]=crane->dX
  f[0] = y[1];
  f[1] = (-(2*c->dr+c->C)*y[1] -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
  //  f[1] = (-2*c->dr -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
  f[2] = y[3];
  f[3] = c->ddX;
  //  if(fabs(t-(int)t)<1e6) fprintf(stderr,"check t=%g f=%g,%g,%g,%g\n ",t,f[0],f[1],f[2],f[3]);
  //  if(fabs(t-(int)(t+0.5))<1e8) fprintf(stderr,"check t=%g y=%g,%g,%g,%g\n ",t,y[0],y[1],y[2],y[3]);
//  if(ddXn>0){
//    fprintf(stderr,"t=%g check?? ddXn=%g=?=%g\n",t,ddXn,c->ddX);
//  }
//  double dd=fabs(t-(int)(t+0.5)); if(t>5 && dd<1e-12){
//    fprintf(stderr,"check t=%g c->ddX=%g, t=%g, it=%d.dd=%g\n ",t,c->ddX,t,(int)(t+0.5),dd);
//  }
  return GSL_SUCCESS;
}
////////////////
//#define square(x) ((x)*(x))
double square(double x){return((x)*(x));}
//double plant(double uu)
double plant1(double uu,CRANE *crane)
{
  //  fprintf(stderr,"check uu=%g, ",uu);
  if(uu>crane->ddXmax) uu=crane->ddXmax;
  else if(uu<-crane->ddXmax) uu=-crane->ddXmax;
  //  fprintf(stderr,"->%g, \n",uu);
  int n;
  for(n=0;n<crane->nh;n++){
    crane->ddX=uu;//crane->ddX=(uu<crane->ddXmax)?uu:crane->ddXmax;
//    if(ddXn>0){
//      fprintf(stderr,"ddXn=%g=?=%g\n",ddXn,crane->ddX);
//    }
    if(crane->dX>=crane->dXmax && crane->ddX>0) crane->ddX=0;
    else if(crane->dX<=-crane->dXmax && crane->ddX<0) crane->ddX=0;
    int status = gsl_odeiv_step_apply (crane->_s, crane->_t, crane->_h,
				       crane->_y, crane->_y_err,
				       crane->_dydt_in,
				       crane->_dydt_out,
				       &crane->_sys);
    if (status != GSL_SUCCESS) {fprintf(stderr,"### failure of plant calc.\n"); break;}

    int i;for(i=0;i<dim_crane;i++) crane->_dydt_in[i]=crane->_dydt_out[i];
    crane->_t+=crane->_h;
    crane->a  =crane->_y[0];
    crane->da =crane->_y[1];
    crane->X  =crane->_y[2];
    crane->dX =crane->_y[3];
    crane->x  =crane->X+crane->r * sin(crane->a); //output y
    crane->y  =crane->r*cos(crane->a);

    crane->dx=(crane->x-crane->x0)/crane->_h;
    crane->dy=(crane->y-crane->y0)/crane->_h;
    crane->ddx=(crane->dx-crane->dx0)/crane->_h;
    crane->ddy=(crane->dy-crane->dy0)/crane->_h;
    crane->T=crane->m*sqrt(square(crane->ddx)+square(crane->ddy-crane->g));
    crane->F=crane->M*crane->ddX-crane->T*sin(crane->a);
    crane->dda=(crane->da-crane->da0)/crane->_h;
    crane->x0=crane->x;
    crane->y0=crane->y;
    crane->dx0=crane->dx;
    crane->dy0=crane->dy;
    crane->da0=crane->da;
  }
#ifndef CRANESUB
  fprintf(fp,"%.7e %.7e %.7e", crane->_t,crane->_y[0],crane->_y[1]);//crane->a,crane->da
  fprintf(fp," %.7e %.7e %.7e" ,crane->X,crane->dX,crane->ddX);
  fprintf(fp," %.7e %.7e %.7e" ,crane->x,crane->dx,crane->ddx);
  fprintf(fp," %.7e %.7e %.7e" ,crane->y,crane->dy,crane->ddy);
  fprintf(fp," %.7e %.7e" ,crane->T,crane->F);
  fprintf(fp,"\n");
#endif //#ifndef CRANESUB
  return(crane->x);//crane->x in [0,AP_y=10?] ==> rt in [0,1]
}
int initialize1(CRANE *crane)
{
  crane->cranefunc=crane1func;
  crane->plant=plant1;
  crane->g=9.8;
  crane->M=100;
  crane->m=20;//10;//20
  crane->r=5;//5;
  crane->C=1;//5;

  crane->dr=crane->ddr=0;
  crane->x0=0;
  crane->y0=crane->r;
  crane->dx0=crane->dy0=crane->da0=0;
  crane->dXmax=1.;//??ddXmax=(dXmax-dX)/dt;//??
  crane->ddXmax=0.2;//??ddXmax=(dXmax-dX)/dt;//??
  crane->Fmax=30.0;//check//  crane->Fmax=30;//check
#ifdef CRANESUB
  crane->h=AP_tS;
  crane->m=_crane_m;
  crane->r=_crane_r;
  crane->dXmax=_crane_dXmax;
  if(_AP_umax>0) AP_u_max=crane->ddXmax=_AP_umax;
  else AP_u_max=crane->ddXmax;
  AP_u_min=-AP_u_max;
  //  starttime=0;
  //  totaltime=50;
  rr=AP_r=_AP_r;//10
  rr_kyoyou=_rr_kyoyou;
  //  p=(double *)malloc(buffsize*sizeof(double));
  C_MODE=11;
  //  iteration=_iteration;
#else
  crane->h=0.01;
#endif
  crane->nh=10;
  crane->_h=crane->h/crane->nh;//0.001
  crane->_dim=4;
  if(crane->_y==NULL){
    crane->_y=(double*)malloc(crane->_dim*sizeof(double));
    crane->_y_err=(double*)malloc(crane->_dim*sizeof(double));
    crane->_dydt_in=(double*)malloc(crane->_dim*sizeof(double));
    crane->_dydt_out=(double*)malloc(crane->_dim*sizeof(double));
  }
  crane->_T= gsl_odeiv_step_rk4;
  crane->_s= gsl_odeiv_step_alloc (crane->_T, crane->_dim);
  crane->_sys= (gsl_odeiv_system){crane1func, NULL, crane->_dim, crane};
  //  crane->_sys= (gsl_odeiv_system){crane1func, NULL, crane->_dim, &crane};
  //  crane->_sys= (gsl_odeiv_system){crane->cranefunc, NULL, crane->_dim, &crane};
  //initialize variables
  crane->_t=crane->_y[0]=crane->_y[1]=crane->_y[2]=crane->_y[3]=0;
  int i;for(i=0;i<crane->_dim;i++) crane->_dydt_in[i]=0;
  GSL_ODEIV_FN_EVAL(&(crane->_sys), crane->_t, crane->_y, crane->_dydt_in);
//  crane->a  =crane->_y[0];
//  crane->da =crane->_y[1];
//  crane->X  =crane->_y[2];
//  crane->dX =crane->_y[3];
//  crane->x  =crane->X+crane->r * sin(crane->a); //output y
//  crane->y  =crane->r*cos(crane->a);
//  crane->dXmax=1.;//??ddXmax=(dXmax-dX)/dt;//??
//  crane->ddXmax=0.2;//??ddXmax=(dXmax-dX)/dt;//??
  return(0);
}//endof initialize1

#ifndef CRANESUB
int main1(int argc,char **argv)
{
  ////////////////////////////////////////////////////
  /// method 1                                     ///
  /// input ddX                                    ///
  /// output a,x,y,F                               ///
  ////////////////////////////////////////////////////
  //  double h = 0.001;//,hh=0.1; h = 0.0001;//,hh=0.1;
  double t0=5;     //stationary
  double t1=5.0+t0;//accelerate(speed up)
  double t2=2.0+t1;//free run
  double t3=5.0+t2;//decelerate(slow down)
  double t4=20+t3;   //free run
  CRANE crane[1];
//  int (*initialize)();
//  if(1==1){
//    initialize=initialize1;
//  }
//  (*initialize)(crane);
  crane->_y=NULL;
  initialize1(crane);
  int M=(t4/crane->h)+1;
  //  ddX=(double*)malloc(sizeof(double)*M);
  double *ddX=(double*)malloc(sizeof(double)*M);
  //  double *F  =(double*)malloc(sizeof(double)*M);
  double Vmax=1;//dX/dt=1[m/s]
  int n;
  double t;
  ////////////////////////////////////////////////////
  ///   set ddX                                    ///
  ////////////////////////////////////////////////////
  int n0 =t0/crane->h+0.5, n4=t4/crane->h+0.5;
  for(t=0;t<t4;t+=crane->h){
    n=t/crane->h;
    if(t<=t0) ddX[n]=0; //zero
    else if(t<=t1) ddX[n]=Vmax/(t1-t0); //accelerate
    else if(t<=t2) ddX[n]=0;            //free move
    else if(t<=t3) ddX[n]=-Vmax/(t3-t2);//decelerate
    else ddX[n]=0;                      //free move
  }
  ////////////////////////////////////////////////////
  ///   Solve by the Runge Kutta Method of GSL     ///
  ////////////////////////////////////////////////////
  //  char *fn="crane1io.dat";
  fp=fopen(fn1,"w");
  fprintf(fp,"#%.7e %d %d %.7e %.7e %.7e %.7e #h,n0,n4,M,n,r\n",crane->h,n0,n4,crane->M,crane->m,crane->r,crane->C);
  //  double dx0=0,dy0=0,da0=0;
  //  for(n=0;n<n4;){
  for(t=0;t<t4;t+=crane->h){
    n=t/crane->h;
//    ddXn=ddX[n];if(n>500 && n<510){
//      fprintf(stderr,"check n=%d,t=%g ddX=%g\n",n,t,ddX[n]);
//    }
    //    (*(crane->plant))(ddX[n],crane);
    plant1(ddX[n],crane);
  }
  fclose(fp);
  fprintf(stdout,"Results are stored in '%s'.\n",fn1);

  gsl_odeiv_step_free (crane->_s);
  return 0;
}//endof main1
#endif//endof #ifndef CRANESUB
/////////////////////////
int crane2func (double t, const double y[], double f[], void *params)
{
  CRANE *c = (CRANE *)params;
  //y[0]=crane->p, y[1]=crane->da, y[2]=crane->X, y[3]=crane->dX
  f[0] = y[1];
  f[1] = (-(2*c->dr+c->C)*y[1] -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
  //  f[1] = (-2*c->dr -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
  f[2] = y[3];
  //  f[3] = c->ddX;
  f[3] = (c->F+c->T*sin(y[0]))/c->M;
  return GSL_SUCCESS;
}//endof crane2func

double plant2(double uu,CRANE *crane)
{
  //  crane->F=uu;
  if(uu>crane->Fmax) uu=crane->Fmax;
  else if(uu<-crane->Fmax) uu=-crane->Fmax;
  int n;
  for(n=0;n<crane->nh;n++){
    crane->F=uu;//
    if(crane->dX>=crane->dXmax && crane->F>0) crane->F=0;
    else if(crane->dX<=-crane->dXmax && crane->F<0) crane->F=0;
    int status = gsl_odeiv_step_apply (crane->_s, crane->_t, crane->_h,
				       crane->_y, crane->_y_err,
				       crane->_dydt_in,
				       crane->_dydt_out,
				       &crane->_sys);
    if (status != GSL_SUCCESS)   break;

    int i;for(i=0;i<dim_crane;i++) crane->_dydt_in[i]=crane->_dydt_out[i];
    crane->_t+=crane->_h;

    crane->a  =crane->_y[0];
    crane->da =crane->_y[1];
    crane->X  =crane->_y[2];
    crane->dX =crane->_y[3];

    crane->x  =crane->X+crane->r * sin(crane->a);
    crane->y  =crane->r*cos(crane->a);
    crane->dx=(crane->x-crane->x0)/crane->_h;
    crane->dy=(crane->y-crane->y0)/crane->_h;
    crane->ddx=(crane->dx-crane->dx0)/crane->_h;
    crane->ddy=(crane->dy-crane->dy0)/crane->_h;
    crane->T=crane->m*sqrt(square(crane->ddx)+square(crane->ddy-crane->g));
    crane->ddX=(crane->dX-crane->dX0)/crane->_h;
    //    crane->F=crane->M*crane->ddX-crane->T*sin(crane->a);
    crane->dda=(crane->da-crane->da0)/crane->_h;
    /*    crane->dda=(-2*crane->dr*crane->da-crane->g*sin(crane->a)-ddX[n]*cos(crane->a))/crane->r; */
    crane->x0=crane->x;
    crane->y0=crane->y;
    crane->dx0=crane->dx;
    crane->dy0=crane->dy;
    crane->da0=crane->da;
    crane->dX0=crane->dX;
  }
#ifndef CRANESUB
  fprintf(fp,"%.7e %.7e %.7e", crane->_t,crane->_y[0],crane->_y[1]);//crane->a,crane->da
  fprintf(fp," %.7e %.7e %.7e" ,crane->X,crane->dX,crane->ddX);
  fprintf(fp," %.7e %.7e %.7e" ,crane->x,crane->dx,crane->ddx);
  fprintf(fp," %.7e %.7e %.7e" ,crane->y,crane->dy,crane->ddy);
  fprintf(fp," %.7e %.7e" ,crane->T,crane->F);
  fprintf(fp,"\n");
#endif
  return(crane->x);
  //  return(crane->X);
}//endof plant2
int initialize2(CRANE *crane)
{
  crane->cranefunc=crane2func;
  crane->plant=plant2;
  crane->Fmax=20.0;//check//  crane->Fmax=30;//check
  //  crane->dXmax=1.0;
  crane->dXmax=1.0;
  crane->ddXmax=0.2;
  crane->g=9.8;
  crane->M=100;
  crane->m=20;//10;//10;//20
  crane->r=5;//5;
  crane->C=1;//5;
  crane->dr=crane->ddr=0;
  crane->x0=0;
  crane->y0=crane->r;
  crane->dX0=crane->dx0=crane->dy0=crane->da0=0;
#ifdef CRANESUB
  //  crane->Fmax=_crane_Fmax;//check//  crane->Fmax=30;//check
  crane->m=_crane_m;
  crane->r=_crane_r;
  crane->h=AP_tS;
  crane->dXmax=_crane_dXmax;
  if(_AP_umax>0) AP_u_max=crane->Fmax=_AP_umax;
  else AP_u_max=crane->Fmax;
  AP_u_min=-AP_u_max;
  //  starttime=0;
  //  totaltime=100;//  totaltime=40;
  //  totaltime=50;//  totaltime=40;
  //  int kmax=totaltime/AP_tS;
  //  _rr=(double*)malloc(kmax*sizeof(double));
  rr=AP_r=_AP_r;//10
  //  int k;for(k=0;k<kmax;k++) _rr[k]=AP_r;//10;
  rr_kyoyou=_rr_kyoyou;
  //  p=(double *)malloc(buffsize*sizeof(double));
  C_MODE=11;
  //  iteration=_iteration;
  //  fprintf(stdout,"?????????iteration=%d\n",_iteration);
#else
  crane->h=0.01;
#endif
  crane->nh=10;
  crane->F=crane->a=crane->da=crane->dda=0;
  crane->T=crane->m*crane->g;
  crane->_h=crane->h/crane->nh;//0.001
  crane->_dim=4;
  {//  if(crane->_y==NULL){
    crane->_y=(double*)malloc(crane->_dim*sizeof(double));
    crane->_y_err=(double*)malloc(crane->_dim*sizeof(double));
    crane->_dydt_in=(double*)malloc(crane->_dim*sizeof(double));
    crane->_dydt_out=(double*)malloc(crane->_dim*sizeof(double));
  }
  crane->_T= gsl_odeiv_step_rk4;
  crane->_s= gsl_odeiv_step_alloc (crane->_T, crane->_dim);
  crane->_sys= (gsl_odeiv_system){crane->cranefunc, NULL, crane->_dim, crane};
  //  crane->_sys= (gsl_odeiv_system){crane->cranefunc, NULL, crane->_dim, &crane};
  crane->_t=crane->_y[0]=crane->_y[1]=crane->_y[2]=crane->_y[3]=0;
  int i;for(i=0;i<crane->_dim;i++) crane->_dydt_in[i]=0;
  GSL_ODEIV_FN_EVAL(&(crane->_sys), crane->_t, crane->_y, crane->_dydt_in);
  return(0);
}
#ifndef CRANESUB
int main2(int argc, char** argv)
{
  //  CRANE crane;
  CRANE crane[1];
  crane->g=9.8;
  crane->M=100;
  crane->m=20;
  crane->r=5;
  crane->dr=crane->ddr=0;
  ////////////////////////////////////////////////////
  /// method 2                                     ///
  /// input F                                     ///
  /// output p,x,y,X,...                          ///
  ////////////////////////////////////////////////////
  double h;//,hh=0.1;
  //result for different h
  /*gnuplot
   set style data lines;n=1;
   n=n+1;plot "crane1io.dat" using 1:n, "crane2io.dat" using 1:n;print "n=",n
   #strange at n=9(ddx),12(ddy),13(T)
  */
  ////////////////////////////////////////////////////
  ///   read F                                    ///
  ////////////////////////////////////////////////////
  //  char *fn1="crane1io.dat";
  FILE *fp1=fopen(fn1,"r");
  int n4;
#define buffsize 512
  char buff[buffsize];
  int n0;
  if(fgets(buff,buffsize,fp1)!=NULL)
    sscanf(buff,"#%lf%d%d%lf%lf%lf%lf",&h,&n0,&n4,&crane->M,&crane->m,&crane->r,&crane->C);
  crane->h=h;
  double t4=n4*h;
  double *F  =(double*)malloc(sizeof(double)*(n4+1));
  //  int M=n4+1;double *F  =(double*)malloc(sizeof(double)*M);
  //  double *ddX=(double*)malloc(sizeof(double)*M);

  int n;
  double t;
  //  double dx0=0,dy0=0,dX0=0,da0=0;
  //  double x0=0,y0=crane->r;

  if(fgets(buff,buffsize,fp1)!=NULL);  //??  fgets(buff,buffsize,fp1);
  for(n=0;n<n4;n++){//
    double _x1,_x2,_x3,_x4,_x5,_x6,_x7,_x8,_x9,_x10,_x11,_x12,_x13,_x14;
    if(fgets(buff,buffsize,fp1)!=NULL) //      fgets(buff,buffsize,fp1);
      sscanf(buff,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
	     &_x1,&_x2,&_x3,&_x4,&_x5,&_x6,&_x7,&_x8,&_x9,&_x10,&_x11,&_x12,&_x13,&_x14);
    F[n]=_x14;
  }
  ////////////////////////////////////////////////////
  ///   Solve by the Runge Kutta Method of GSL     ///
  ////////////////////////////////////////////////////
  crane->_y=NULL;
  initialize2(crane);
  fp=fopen(fn2,"w");
  for(t=0;t<t4;t+=crane->h){
    n=t/crane->h;
    plant2(F[n],crane);
    //    (*(crane->plant))(F[n],crane);
  }
  fclose(fp);
  fprintf(stdout,"Results are stored in '%s'.\n",fn2);

  gsl_odeiv_step_free (crane->_s);
  return 0;
}//main2
#endif //#ifndef CRANESUB

#ifndef CRANESUB
int main(int argc,char **argv)
{
  if(argc<2){
    fprintf(stderr,"Usage:%s [1|2]\n",argv[0]);
    return(1);
  }
  if(argv[1][0]!='2') main1(argc,argv);
  else main2(argc,argv);
  return(0);
}
#endif //#ifndef CRANESUB


text:
/*
copied from http://www.gnu.org/software/gsl/manual/html_node/Basis-Splines.html
*/

#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_sf_erf.h>

void geterf_(double *xanderfx) {
  xanderfx[1] = gsl_sf_erf(xanderfx[0]);
}

//spline locations held fixed in Mpc^-1; CMB basically fixed P(k) in these units
void dopksmoothbspline_(double *kvals, double *lnpklinear, double *lnpksmooth, int *npts)  {

	double kmaxsuppress = 0.01*0.7;
	size_t n, ncoeffs, nbreak;
	gsl_bspline_workspace *bw;
	gsl_vector *B;
	gsl_vector *c, *w, *x, *y;
	gsl_matrix *X, *cov;
	gsl_multifit_linear_workspace *mw;
	double deltak,lastk;
	int i,j,countkeep;

	nbreak = 9;
	gsl_vector *mybreaks = gsl_vector_alloc(nbreak);
	gsl_vector_set(mybreaks,0,(0.001*0.7));
	gsl_vector_set(mybreaks,1,(0.025*0.7));
	gsl_vector_set(mybreaks,2,(0.075*0.7));
	gsl_vector_set(mybreaks,3,(0.125*0.7));
	gsl_vector_set(mybreaks,4,(0.175*0.7));
	gsl_vector_set(mybreaks,5,(0.225*0.7));
	gsl_vector_set(mybreaks,6,(0.275*0.7));
	gsl_vector_set(mybreaks,7,(0.325*0.7));
	gsl_vector_set(mybreaks,8,(0.375*0.7));

	countkeep = 0;
	for(i=0;i<(*npts);i++)  {
		if((kvals[i]) >= gsl_vector_get(mybreaks,0) && (kvals[i]) <= gsl_vector_get(mybreaks,nbreak-1)) {
			countkeep += 1;
			}
		}
	n = countkeep;
	ncoeffs = nbreak + 2;

	/* allocate a cubic bspline workspace (k = 4) */
	bw = gsl_bspline_alloc(4, nbreak);
	B = gsl_vector_alloc(ncoeffs);
	x = gsl_vector_alloc(n);
	y = gsl_vector_alloc(n);
	X = gsl_matrix_alloc(n, ncoeffs);
	c = gsl_vector_alloc(ncoeffs);
	w = gsl_vector_alloc(n);
	cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
	mw = gsl_multifit_linear_alloc(n, ncoeffs);
	i=0;
	for(j=0;j<(*npts);j++)  {
		if((kvals[j]) >= gsl_vector_get(mybreaks,0) && (kvals[j]) <= gsl_vector_get(mybreaks,nbreak-1)) {
			gsl_vector_set(x,i,(kvals[j]));
			gsl_vector_set(y,i,exp(lnpklinear[j])*pow(kvals[j],1.5));
			if(j>0)  {
				deltak = kvals[j] - kvals[j-1];
				}
			else {
				deltak = kvals[0];
				if(kvals[1] - kvals[0] < deltak)  {
					deltak = kvals[1]-kvals[0];
					}
				}
			gsl_vector_set(w,i,deltak);
			i+=1;
			}
		}
	gsl_bspline_knots(mybreaks,bw);
	for(i=0;i<n;i++)  {
		double xi = gsl_vector_get(x,i);
		gsl_bspline_eval(xi,B,bw);
		for(j=0;j<ncoeffs;j++)  {
			double Bj = gsl_vector_get(B,j);
			gsl_matrix_set(X,i,j,Bj);
			}
		}
	//do fit
	double yi,yierr,chisq;
	gsl_multifit_wlinear(X,w,y,c,cov,&chisq,mw);
	i = 0;
	for(j=0;j<(*npts);j++)  {
		if((kvals[j]) >= gsl_vector_get(mybreaks,0) && (kvals[j]) <= gsl_vector_get(mybreaks,nbreak-1)) {
			gsl_bspline_eval(gsl_vector_get(x,i),B,bw);
			gsl_multifit_linear_est(B,c,cov,&yi,&yierr);
			lnpksmooth[j] = log(yi*pow(kvals[j],-1.5));
			i += 1;
			}
		else {
			lnpksmooth[j] = lnpklinear[j];
			}
		//spline is wacky at small k -- suppress difference at k < 0.01
		if(kvals[j] < kmaxsuppress)  {
			lnpksmooth[j] = lnpklinear[j];
			}
		}
	assert(i==n);
	gsl_bspline_free(bw);
	gsl_vector_free(B);
	gsl_vector_free(x);
	gsl_vector_free(y);
	gsl_vector_free(mybreaks);
	gsl_matrix_free(X);
	gsl_vector_free(c);
	gsl_vector_free(w);
	gsl_matrix_free(cov);
	gsl_multifit_linear_free(mw);
	}
/*

FILE *open_file_read(char *filename) {
	FILE *ifp;

	if(!(ifp=fopen(filename,"r")))  {
		fprintf(stderr,"Can't open %s\n",filename);
		return NULL;
		}
	else {
		return ifp;
		}
	}

int get_file_length(FILE *ifp, int *headercount, int *linecount) {
	int i;
	char line[256];

	*headercount = 0;
	*linecount = 0;

    while(fgets(line,256,ifp))  {
        if(!feof(ifp))  {
            (*linecount)++;
            }
		else {
			break;
			}
		if(strstr(line,"#"))  {
			(*headercount)++;
			}
        }
	(void) rewind(ifp);
	}

int main()  {
	double *lnpksmooth, *lnpklinear, *kvals;
	FILE *ifp;
	int headercount, linecount,i,j;
	double kval,pnw,plin;
	ifp = open_file_read("/Users/breid/cosmomcmnu/test1.cut");
	get_file_length(ifp,&headercount,&linecount);
	lnpksmooth = (double *) malloc(sizeof(double)*linecount);
	lnpklinear = (double *) malloc(sizeof(double)*linecount);
	kvals = (double *) malloc(sizeof(double)*linecount);
	for(i=0;i<linecount;i++)  {
		fscanf(ifp,"%lf %lf %lf\n",&kval,&pnw,&plin);
		lnpklinear[i] = plin;
		kvals[i] = kval*0.701;
		//kvals[i] = kval;
		}
	dopksmoothbspline_(kvals,lnpklinear,lnpksmooth,&linecount);
	for(i=0;i<linecount;i++)  {
		printf("%f %f %f\n",kvals[i],lnpklinear[i],lnpksmooth[i]);
		}
	return 0;
	}
*/

max_stars_repo_name: sjoudaki/CosmoJBD

max_issues_repo_name: sjoudaki/CosmoJBD

max_forks_repo_name: sjoudaki/CosmoJBD