text:
//
//  trmm.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 1/4/12.
//  Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//

#ifndef _trmm_h
#define _trmm_h

/// @file trmm.h Performs multiplication of a triangular matrix with a rectangular matrix.

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

namespace LATL
{
   /// @brief Performs multiplication of a real triangular matrix with a real rectangular matrix.
   ///
   /// For a real upper or lower triangular matrix A, real rectangular matrix B and real scalar alpha,
   ///
   ///        B := alpha*A*B  or  B := alpha*A'*B  or  B := alpha*B*A  or  B := alpha*B*A'
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param side Specifies whether the matrix A appears on the left or right side as follows:
   ///
   ///        if side = 'L' or 'l' then B := alpha*A*B or alpha*A'*B
   ///        if side = 'R' or 'r' then C := alpha*B*A or alpha*B*A'
   /// @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 trans Specifies wheather the transpose of A is to be used or not:
   ///
   ///        if trans = 'N' or 'n' then B := alpha*A*B or alpha*B*A
   ///        if trans = 'T' or 't' then B := alpha*A'*B or alpha*B*A'
   ///        if trans = 'C' or 'c' then B := alpha*A'*B or alpha*B*A'
   /// @param diag specifies whether or not A is unit triangular as follows:
   ///
   ///        if diag = 'U' or 'u' then A is assumed to be unit triangular
   ///        if diag = 'N' or 'n' then A is not assumed to be unit triangular.
   /// If A is unit triangular, the diagonal elements are assumed to be unity and are not referenced.
   /// @param m Specifies the number of rows of the matrix B.  m>=0
   /// @param n Specifies the number of columns of the matrix B.  n>=0
   /// @param alpha Real scalar.
   /// @param A Pointer to real triangular matrix A.  The order of A is n if side = 'L' or 'l';
   /// the order of A is m is side = 'R' or 'r'.
   /// If uplo = 'U' or 'u, A is upper triangular and the lower triangular part is not referenced.
   /// If uplo = 'L' or 'l A is lower triangular and the upper triangular part is not referenced.
   /// @param ldA Column length of the matrix A.  If side='L' or 'l' then ldA>=n; if side='R' or 'r' then lda>=m.
   /// @param B Pointer to real m-by-n matrix B.
   /// @param ldB Column length of the matrix B.  ldB>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int TRMM(char side, char uplo, char trans, char diag, int_t m, int_t n, real_t alpha, real_t *A, int_t ldA, real_t *B, int_t ldB)
   {
      using std::toupper;

      const real_t zero(0.0);
      int_t i,j,k;
      real_t *a,*b,*bt;
      real_t t;

      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      bool nounit=(diag=='N')?1:0;

      if(alpha==zero)
      {
         b=B;
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               b[i]=zero;
            b+=ldB;
         }
      }
      else if(side=='L')
      {
         if(trans=='N')
         {
            if(uplo=='U')
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A;
                  for(k=0;k<m;k++)
                  {
                     t=alpha*b[k];
                     for(i=0;i<k;i++)
                        b[i]+=t*a[i];
                     if(nounit)
                        t*=a[k];
                     b[k]=t;
                     a+=ldA;
                  }
                  b+=ldB;
               }
            }
            else
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A+m*ldA;
                  for(k=m-1;k>=0;k--)
                  {
                     a-=ldA;
                     t=alpha*b[k];
                     b[k]=t;
                     if(nounit)
                        b[k]*=a[k];
                     for(i=k+1;i<m;i++)
                        b[i]+=t*a[i];
                  }
                  b+=ldB;
               }
            }
         }
         else
         {
            if(uplo=='U')
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A+m*ldA;
                  for(i=m-1;i>=0;i--)
                  {
                     a-=ldA;
                     t=b[i];
                     if(nounit)
                        t*=a[i];
                     for(k=0;k<i;k++)
                        t+=a[k]*b[k];
                     b[i]=alpha*t;
                  }
                  b+=ldB;
               }
            }
            else
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A;
                  for(i=0;i<m;i++)
                  {
                     t=b[i];
                     if(nounit)
                        t*=a[i];
                     for(k=i+1;k<m;k++)
                        t+=a[k]*b[k];
                     b[i]=alpha*t;
                     a+=ldA;
                  }
                  b+=ldB;
               }
            }
         }
      }
      else
      {
         if(trans=='N')
         {
            if(uplo=='U')
            {
               a=A+n*ldA;
               b=B+n*ldB;
               for(j=n-1;j>=0;j--)
               {
                  a-=ldA;
                  b-=ldB;
                  if(nounit)
                     t=alpha*a[j];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  bt=B;
                  for(k=0;k<j;k++)
                  {
                     t=alpha*a[k];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     bt+=ldB;
                  }
               }
            }
            else
            {
               a=A;
               b=B;
               for(j=0;j<n;j++)
               {
                  if(nounit)
                     t=alpha*a[j];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  bt=b+ldB;
                  for(k=j+1;k<n;k++)
                  {
                     t=alpha*a[k];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     bt+=ldB;
                  }
                  a+=ldA;
                  b+=ldB;
               }
            }
         }
         else
         {
            if(uplo=='U')
            {
               a=A;
               bt=B;
               for(k=0;k<n;k++)
               {
                  b=B;
                  for(j=0;j<k;j++)
                  {
                     t=alpha*a[j];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     b+=ldB;
                  }
                  if(nounit)
                     t=alpha*a[k];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  a+=ldA;
                  bt+=ldB;
               }
            }
            else
            {
               a=A+n*ldA;
               bt=B+n*ldB;
               for(k=n-1;k>=0;k--)
               {
                  a-=ldA;
                  bt-=ldB;
                  b=B+(k+1)*ldB;
                  for(j=k+1;j<n;j++)
                  {
                     t=alpha*a[j];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     b+=ldB;
                  }
                  if(nounit)
                     t=alpha*a[k];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     bt[i]*=t;
               }
            }
         }
      }
      return 0;
   }

   /// @brief Performs multiplication of a complex triangular matrix with a complex rectangular matrix.
   ///
   /// For a complex upper or lower triangular matrix A, complex rectangular matrix B and complex scalar alpha,
   ///
   ///        B:=alpha*A*B or B:=alpha*A'*B or B:=alpha*A.'*B or B:=alpha*B*A or B:=alpha*B*A' or B:=alpha*B*A.'
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param side Specifies whether the matrix A appears on the left or right side as follows:
   ///
   ///        if side = 'L' or 'l' then B := alpha*A*B or alpha*A'*B or alpha*A.'*B
   ///        if side = 'R' or 'r' then C := alpha*B*A or alpha*B*A' or alpha*B*A.'
   /// @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 trans Specifies wheather the transpose or conjugate transpose is applied to A as follows:
   ///
   ///        if trans = 'N' or 'n' then B := alpha*A*B or alpha*B*A
   ///        if trans = 'T' or 't' then B := alpha*A.'*B or alpha*B*A.'
   ///        if trans = 'C' or 'c' then B := alpha*A'*B or alpha*B*A'
   /// @param diag specifies whether or not A is unit triangular as follows:
   ///
   ///        if diag = 'U' or 'u' then A is assumed to be unit triangular
   ///        if diag = 'N' or 'n' then A is not assumed to be unit triangular.
   /// If A is unit triangular, the diagonal elements are assumed to be unity and are not referenced.
   /// @param m Specifies the number of rows of the matrix B.  m>=0
   /// @param n Specifies the number of columns of the matrix B.  n>=0
   /// @param alpha Complex scalar.
   /// @param A Pointer to complex triangular matrix A.  The order of A is n if side = 'L' or 'l';
   /// the order of A is m is side = 'R' or 'r'.
   /// If uplo = 'U' or 'u, A is upper triangular and the lower triangular part is not referenced.
   /// If uplo = 'L' or 'l A is lower triangular and the upper triangular part is not referenced.
   /// @param ldA Column length of the matrix A.  If side='L' or 'l' then ldA>=n; if side='R' or 'r' then lda>=m.
   /// @param B Pointer to complex m-by-n matrix B.
   /// @param ldB Column length of the matrix B.  ldB>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int TRMM(char side, char uplo, char trans, char diag, int_t m, int_t n, complex<real_t> alpha, complex<real_t> *A, int_t ldA, complex<real_t> *B, int_t ldB)
   {
      using std::conj;
      using std::toupper;

      const complex<real_t> zero(0.0,0.0);
      int_t i,j,k;
      complex<real_t> *a,*b,*bt;
      complex<real_t> t;

      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      bool nounit=(diag=='N')?1:0;

      if(alpha==zero)
      {
         b=B;
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               b[i]=zero;
            b+=ldB;
         }
      }
      else if(side=='L')
      {
         if(trans=='N')
         {
            if(uplo=='U')
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A;
                  for(k=0;k<m;k++)
                  {
                     t=alpha*b[k];
                     for(i=0;i<k;i++)
                        b[i]+=t*a[i];
                     if(nounit)
                        t*=a[k];
                     b[k]=t;
                     a+=ldA;
                  }
                  b+=ldB;
               }
            }
            else
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A+m*ldA;
                  for(k=m-1;k>=0;k--)
                  {
                     a-=ldA;
                     t=alpha*b[k];
                     b[k]=t;
                     if(nounit)
                        b[k]*=a[k];
                     for(i=k+1;i<m;i++)
                        b[i]+=t*a[i];
                  }
                  b+=ldB;
               }
            }
         }
         else if(trans=='T')
         {
            if(uplo=='U')
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A+m*ldA;
                  for(i=m-1;i>=0;i--)
                  {
                     a-=ldA;
                     t=b[i];
                     if(nounit)
                        t*=a[i];
                     for(k=0;k<i;k++)
                        t+=a[k]*b[k];
                     b[i]=alpha*t;
                  }
                  b+=ldB;
               }
            }
            else
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A;
                  for(i=0;i<m;i++)
                  {
                     t=b[i];
                     if(nounit)
                        t*=a[i];
                     for(k=i+1;k<m;k++)
                        t+=a[k]*b[k];
                     b[i]=alpha*t;
                     a+=ldA;
                  }
                  b+=ldB;
               }
            }
         }
         else if(trans=='C')
         {
            if(uplo=='U')
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A+m*ldA;
                  for(i=m-1;i>=0;i--)
                  {
                     a-=ldA;
                     t=b[i];
                     if(nounit)
                        t*=conj(a[i]);
                     for(k=0;k<i;k++)
                        t+=conj(a[k])*b[k];
                     b[i]=alpha*t;
                  }
                  b+=ldB;
               }
            }
            else
            {
               b=B;
               for(j=0;j<n;j++)
               {
                  a=A;
                  for(i=0;i<m;i++)
                  {
                     t=b[i];
                     if(nounit)
                        t*=conj(a[i]);
                     for(k=i+1;k<m;k++)
                        t+=conj(a[k])*b[k];
                     b[i]=alpha*t;
                     a+=ldA;
                  }
                  b+=ldB;
               }
            }
         }
      }
      else
      {
         if(trans=='N')
         {
            if(uplo=='U')
            {
               a=A+n*ldA;
               b=B+n*ldB;
               for(j=n-1;j>=0;j--)
               {
                  a-=ldA;
                  b-=ldB;
                  if(nounit)
                     t=alpha*a[j];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  bt=B;
                  for(k=0;k<j;k++)
                  {
                     t=alpha*a[k];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     bt+=ldB;
                  }
               }
            }
            else
            {
               a=A;
               b=B;
               for(j=0;j<n;j++)
               {
                  if(nounit)
                     t=alpha*a[j];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  bt=b+ldB;
                  for(k=j+1;k<n;k++)
                  {
                     t=alpha*a[k];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     bt+=ldB;
                  }
                  a+=ldA;
                  b+=ldB;
               }
            }
         }
         else if(trans=='T')
         {
            if(uplo=='U')
            {
               a=A;
               bt=B;
               for(k=0;k<n;k++)
               {
                  b=B;
                  for(j=0;j<k;j++)
                  {
                     t=alpha*a[j];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     b+=ldB;
                  }
                  if(nounit)
                     t=alpha*a[k];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  a+=ldA;
                  bt+=ldB;
               }
            }
            else
            {
               a=A+n*ldA;
               bt=B+n*ldB;
               for(k=n-1;k>=0;k--)
               {
                  a-=ldA;
                  bt-=ldB;
                  b=B+(k+1)*ldB;
                  for(j=k+1;j<n;j++)
                  {
                     t=alpha*a[j];
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     b+=ldB;
                  }
                  if(nounit)
                     t=alpha*a[k];
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     bt[i]*=t;
               }
            }
         }
         else if(trans=='C')
         {
            if(uplo=='U')
            {
               a=A;
               bt=B;
               for(k=0;k<n;k++)
               {
                  b=B;
                  for(j=0;j<k;j++)
                  {
                     t=alpha*conj(a[j]);
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     b+=ldB;
                  }
                  if(nounit)
                     t=alpha*conj(a[k]);
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     b[i]*=t;
                  a+=ldA;
                  bt+=ldB;
               }
            }
            else
            {
               a=A+n*ldA;
               bt=B+n*ldB;
               for(k=n-1;k>=0;k--)
               {
                  a-=ldA;
                  bt-=ldB;
                  b=B+(k+1)*ldB;
                  for(j=k+1;j<n;j++)
                  {
                     t=alpha*conj(a[j]);
                     for(i=0;i<m;i++)
                        b[i]+=t*bt[i];
                     b+=ldB;
                  }
                  if(nounit)
                     t=alpha*conj(a[k]);
                  else
                     t=alpha;
                  for(i=0;i<m;i++)
                     bt[i]*=t;
               }
            }
         }
      }
      return 0;
   }

#ifdef __latl_cblas
#include <cblas.h>

   template <> int TRMM<float>(char side, char uplo, char trans, char diag, int_t m, int_t n, float alpha, float *A, int_t ldA, float *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:CblasTrans;

      cblas_strmm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,alpha,A,ldA,B,ldB);

      return 0;
   }

   template <> int TRMM<double>(char side, char uplo, char trans, char diag, int_t m, int_t n, double alpha, double *A, int_t ldA, double *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:CblasTrans;

      cblas_dtrmm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,alpha,A,ldA,B,ldB);

      return 0;
   }
   template <> int TRMM<float>(char side, char uplo, char trans, char diag, int_t m, int_t n, complex<float> alpha, complex<float> *A, int_t ldA, complex<float> *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_ctrmm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,&alpha,A,ldA,B,ldB);

      return 0;
   }
   template <> int TRMM<double>(char side, char uplo, char trans, char diag, int_t m, int_t n, complex<double> alpha, complex<double> *A, int_t ldA, complex<double> *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_ztrmm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,&alpha,A,ldA,B,ldB);

      return 0;
   }
#endif

}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_bessel.h"

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC besi0 */

/* chebyshev expansions

 series for bi0        on the interval  0.	    to  9.00000d+00
					with weighted error   2.46e-18
					 log weighted error  17.61
			       significant figures required  17.90
				    decimal places required  18.15

 series for ai0        on the interval  1.25000d-01 to  3.33333d-01
					with weighted error   7.87e-17
					 log weighted error  16.10
			       significant figures required  14.69
				    decimal places required  16.76

 series for ai02       on the interval  0.	    to  1.25000d-01
					with weighted error   3.79e-17
					 log weighted error  16.42
			       significant figures required  14.86
				    decimal places required  17.09
*/

static double bi0_data[12] = {
  -.07660547252839144951,
  1.92733795399380827000,
   .22826445869203013390,
   .01304891466707290428,
   .00043442709008164874,
   .00000942265768600193,
   .00000014340062895106,
   .00000000161384906966,
   .00000000001396650044,
   .00000000000009579451,
   .00000000000000053339,
   .00000000000000000245
};
static cheb_series bi0_cs = {
  bi0_data,
  11,
  -1, 1,
  11
};

static double ai0_data[21] = {
   .07575994494023796,
   .00759138081082334,
   .00041531313389237,
   .00001070076463439,
  -.00000790117997921,
  -.00000078261435014,
   .00000027838499429,
   .00000000825247260,
  -.00000001204463945,
   .00000000155964859,
   .00000000022925563,
  -.00000000011916228,
   .00000000001757854,
   .00000000000112822,
  -.00000000000114684,
   .00000000000027155,
  -.00000000000002415,
  -.00000000000000608,
   .00000000000000314,
  -.00000000000000071,
   .00000000000000007
};
static cheb_series ai0_cs = {
  ai0_data,
  20,
  -1, 1,
  13
};

static double ai02_data[22] = {
   .05449041101410882,
   .00336911647825569,
   .00006889758346918,
   .00000289137052082,
   .00000020489185893,
   .00000002266668991,
   .00000000339623203,
   .00000000049406022,
   .00000000001188914,
  -.00000000003149915,
  -.00000000001321580,
  -.00000000000179419,
   .00000000000071801,
   .00000000000038529,
   .00000000000001539,
  -.00000000000004151,
  -.00000000000000954,
   .00000000000000382,
   .00000000000000176,
  -.00000000000000034,
  -.00000000000000027,
   .00000000000000003
};
static cheb_series ai02_cs = {
  ai02_data,
  21,
  -1, 1,
  11
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result)
{
  double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - y;
    result->err = 0.5*y*y;
    return GSL_SUCCESS;
  }
  else if(y <= 3.0) {
    const double ey = exp(-y);
    gsl_sf_result c;
    cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
    result->val = ey * (2.75 + c.val);
    result->err = GSL_DBL_EPSILON * fabs(result->val) + ey * c.err;
    return GSL_SUCCESS;
  }
  else if(y <= 8.0) {
    const double sy = sqrt(y);
    gsl_sf_result c;
    cheb_eval_e(&ai0_cs, (48.0/y-11.0)/5.0, &c);
    result->val  = (0.375 + c.val) / sy;
    result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
    result->err += c.err / sy;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sy = sqrt(y);
    gsl_sf_result c;
    cheb_eval_e(&ai02_cs, 16.0/y-1.0, &c);
    result->val = (0.375 + c.val) / sy;
    result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
    result->err += c.err / sy;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result)
{
  double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0;
    result->err = 0.5*y*y;
    return GSL_SUCCESS;
  }
  else if(y <= 3.0) {
    gsl_sf_result c;
    cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
    result->val  = 2.75 + c.val;
    result->err  = GSL_DBL_EPSILON * (2.75 + fabs(c.val));
    result->err += c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(y < GSL_LOG_DBL_MAX - 1.0) {
    const double ey = exp(y);
    gsl_sf_result b_scaled;
    gsl_sf_bessel_I0_scaled_e(x, &b_scaled);
    result->val  = ey * b_scaled.val;
    result->err  = ey * b_scaled.err + y*GSL_DBL_EPSILON*fabs(result->val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_I0_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_I0_scaled_e(x, &result); )
}

double gsl_sf_bessel_I0(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_I0_e(x, &result); )
}

ext: c


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <lapacke.h>
#define ABS(a)     (((a) < 0) ? -(a) : (a))

int main()
{
  /*
   *	hbar^2/2m = 1
  */

  static double pi = M_PI;
  int n, npw;
  double v0, a, b, ecut, k;
  double *g, *e, *h, *wrk;
  int i, j, ij, m, lwork, info, p;
  char *V = "V";
  char *U = "U";
  FILE *out;

  /*  Input data
      Potential: V(x)=-V_0 for |x|<b/2, V(x)=0 for |x|>b/2
      Periodicity:    V(x+a)=V(x)
  */
  fprintf(stdout, "Parameters for potential well: V_0, a, b > ");
  scanf("%lf %lf %lf",&v0, &a, &b);
  if ( v0 <= 0 ||  a <= 0 || b <= 0 || a <= b )  {
    fprintf(stderr, "wrong input parameters\n");
    exit(1);
  }
  fprintf (stdout,"   V_0=%f, a0=%f, b=%f\n",v0,a,b);
  /*
    Plane waves basis set: G_n=n*2pi/a, \hbar^2/2m*G^2 < Ecut
  */
  fprintf(stdout, "Cutoff for plane waves: ecut > ");
  scanf("%lf",&ecut);
  if ( ecut <= 0) {
    fprintf(stderr, "wrong input parameter\n");
    exit(2);
  }
  /*
    Number of plane waves
  */
  npw = (int) ( sqrt ( ecut/pow( 2.0*pi/a, 2) ) + 0.5 );
  npw = 2*npw+1;
  fprintf (stdout,"   ecut = %f,  # of PWs=%d\n",ecut, npw);
  /*
    Assign values of  G_n: n=0,+1,-1,+2,-2, etc
  */
  n=100;
  g  = (double *) malloc ( npw * sizeof (double) );
  e  = (double *) malloc ( npw * sizeof (double) );
  h  = (double *) malloc ( n*n * sizeof (double) );
  wrk= (double *) malloc (3*npw* sizeof (double) );

  g[0] = 0.0;
  for ( i = 1; i < npw; i+=2 ) {
    g[i  ] = (i+1)*pi/a;
    g[i+1] =-(i+1)*pi/a;
  }
  for ( ij = 0; ij < n*n; ++ij ) {
    h[ij] = 0.0;
  }
  /*
    Loop on k-vectors: k runs from -pi/a to pi/a
  */
  out = fopen("bands.out", "w");
  for ( m =-n; m <= n; m++ ) {
    k = m*pi/n/a;
    /*
      Assign values of the matrix elements of the hamiltonian
      on the plane wave basis
    */
    ij = 0;
    for ( i = 0; i < npw; ++i ) {
      for ( j = 0; j < npw; ++j ) {
      /* NOTA BENE: the matrix h is a vector in the calling program,
         while dsyev expects a (pointer to a) fortran matrix.
         A fortran matrix is "simulated" in the following way:
         if h[ij] is the C array and  h(i,j) is a N*M Fortran matrix,
         h[ij] == h(i,j), where  ij = (j-1)*N + (i-1)  (i=1,N, j=1,M)
      */

	if ( i == j ) {
	  h[ij++] = (k+g[i])*(k+g[i]) - v0/a*b;
	} else {
	  h[ij++] = -v0/a * sin( (g[j]-g[i])*b/2.0 ) / (g[j]-g[i])*2.0;
	}
      }
    }
    /*
      Solution [expansion coefficients are stored into h(i,j)
      j=basis function index, i= eigenvalue index]
      (beware fortran-C reversed index convention!)
    */
    lwork = 3*npw;
    dsyev_ ( V, U, &npw, h, &npw, e, wrk, &lwork, &info );

    if ( info != 0) {
      fprintf(stderr, "H-matrix diagonalization failed\n");
      exit(3);
    }

    /*printf("%f ",k);
    for(p=0; p<npw; p++) printf("%f ",e[p]);printf("\n");*/

    fprintf(out,"%12.6f ",k);
    for(p=0; p<npw; p++) fprintf(out,"%10.6f ",e[p]);
    fprintf(out,"\n");

  }
  fclose(out);

  return 0;

}

max_stars_repo_name: algebrato/Struttura

max_issues_repo_name: algebrato/Struttura

max_forks_repo_name: algebrato/Struttura


text:
/*
  C code for the adiabatic approximation
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_integration.h>
#ifdef _OPENMP
#include <omp.h>
#endif
#define CHUNKSIZE 10
//Potentials
#include <galpy_potentials.h>
#include <actionAngle.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/*
  Structure Declarations
*/
struct JRAdiabaticArg{
  double ER;
  double Lz22;
  int nargs;
  struct potentialArg * actionAngleArgs;
};
struct JzAdiabaticArg{
  double Ez;
  double R;
  int nargs;
  struct potentialArg * actionAngleArgs;
};
/*
  Function Declarations
*/
void actionAngleAdiabatic_actions(int,double *,double *,double *,double *,
				 double *,int,int *,double *,double,
				 double *,double *,int *);
void calcJRAdiabatic(int,double *,double *,double *,double *,double *,
		     int,struct potentialArg *,int);
void calcJzAdiabatic(int,double *,double *,double *,double *,int,
		     struct potentialArg *,int);
void calcRapRperi(int,double *,double *,double *,double *,double *,
		  int,struct potentialArg *);
void calcZmax(int,double *,double *,double *,double *,int,
	      struct potentialArg *);
double JRAdiabaticIntegrandSquared(double,void *);
double JRAdiabaticIntegrand(double,void *);
double JzAdiabaticIntegrandSquared(double,void *);
double JzAdiabaticIntegrand(double,void *);
double evaluateVerticalPotentials(double, double,int, struct potentialArg *);
/*
  Actual functions, inlines first
*/
inline void calcEREzL(int ndata,
		      double *R,
		      double *vR,
		      double *vT,
		      double *z,
		      double *vz,
		      double *ER,
		      double *Ez,
		      double *Lz,
		      int nargs,
		      struct potentialArg * actionAngleArgs){
  int ii;
  UNUSED int chunk= CHUNKSIZE;
#pragma omp parallel for schedule(static,chunk) private(ii)
  for (ii=0; ii < ndata; ii++){
    *(ER+ii)= evaluatePotentials(*(R+ii),0.,
				 nargs,actionAngleArgs)
      + 0.5 * *(vR+ii) * *(vR+ii)
      + 0.5 * *(vT+ii) * *(vT+ii);
    *(Ez+ii)= evaluateVerticalPotentials(*(R+ii),*(z+ii),
					 nargs,actionAngleArgs)
      + 0.5 * *(vz+ii) * *(vz+ii);
    *(Lz+ii)= *(R+ii) * *(vT+ii);
  }
}
/*
  MAIN FUNCTIONS
 */
void actionAngleAdiabatic_actions(int ndata,
				  double *R,
				  double *vR,
				  double *vT,
				  double *z,
				  double *vz,
				  int npot,
				  int * pot_type,
				  double * pot_args,
				  double gamma,
				  double *jr,
				  double *jz,
				  int * err){
  int ii;
  //Set up the potentials
  struct potentialArg * actionAngleArgs= (struct potentialArg *) malloc ( npot * sizeof (struct potentialArg) );
  parse_actionAngleArgs(npot,actionAngleArgs,pot_type,pot_args);
  //ER, Ez, Lz
  double *ER= (double *) malloc ( ndata * sizeof(double) );
  double *Ez= (double *) malloc ( ndata * sizeof(double) );
  double *Lz= (double *) malloc ( ndata * sizeof(double) );
  calcEREzL(ndata,R,vR,vT,z,vz,ER,Ez,Lz,npot,actionAngleArgs);
  //Calculate peri and apocenters
  double *rperi= (double *) malloc ( ndata * sizeof(double) );
  double *rap= (double *) malloc ( ndata * sizeof(double) );
  double *zmax= (double *) malloc ( ndata * sizeof(double) );
  calcZmax(ndata,zmax,z,R,Ez,npot,actionAngleArgs);
  calcJzAdiabatic(ndata,jz,zmax,R,Ez,npot,actionAngleArgs,10);
  //Adjust planar effective potential for gamma
  UNUSED int chunk= CHUNKSIZE;
#pragma omp parallel for schedule(static,chunk) private(ii)
  for (ii=0; ii < ndata; ii++){
    *(Lz+ii)= fabs( *(Lz+ii) ) + gamma * *(jz+ii);
    *(ER+ii)+= 0.5 * *(Lz+ii) * *(Lz+ii) / *(R+ii) / *(R+ii)
      - 0.5 * *(vT+ii) * *(vT+ii);
  }
  calcRapRperi(ndata,rperi,rap,R,ER,Lz,npot,actionAngleArgs);
  calcJRAdiabatic(ndata,jr,rperi,rap,ER,Lz,npot,actionAngleArgs,10);
  for (ii=0; ii < npot; ii++) {
    if ( (actionAngleArgs+ii)->i2d )
      interp_2d_free((actionAngleArgs+ii)->i2d) ;
    if ((actionAngleArgs+ii)->accx )
      gsl_interp_accel_free ((actionAngleArgs+ii)->accx);
    if ((actionAngleArgs+ii)->accy )
      gsl_interp_accel_free ((actionAngleArgs+ii)->accy);
    free((actionAngleArgs+ii)->args);
  }
  free(actionAngleArgs);
  free(ER);
  free(Ez);
  free(Lz);
  free(rperi);
  free(rap);
  free(zmax);
}
void calcJRAdiabatic(int ndata,
		     double * jr,
		     double * rperi,
		     double * rap,
		     double * ER,
		     double * Lz,
		     int nargs,
		     struct potentialArg * actionAngleArgs,
		     int order){
  int ii, tid, nthreads;
#ifdef _OPENMP
  nthreads = omp_get_max_threads();
#else
  nthreads = 1;
#endif
  gsl_function * JRInt= (gsl_function *) malloc ( nthreads * sizeof(gsl_function) );
  struct JRAdiabaticArg * params= (struct JRAdiabaticArg *) malloc ( nthreads * sizeof (struct JRAdiabaticArg) );
  for (tid=0; tid < nthreads; tid++){
    (params+tid)->nargs= nargs;
    (params+tid)->actionAngleArgs= actionAngleArgs;
  }
  //Setup integrator
  gsl_integration_glfixed_table * T= gsl_integration_glfixed_table_alloc (order);
  UNUSED int chunk= CHUNKSIZE;
#pragma omp parallel for schedule(static,chunk)				\
  private(tid,ii)							\
  shared(jr,rperi,rap,JRInt,params,T,ER,Lz)
  for (ii=0; ii < ndata; ii++){
#ifdef _OPENMP
    tid= omp_get_thread_num();
#else
    tid = 0;
#endif
    if ( *(rperi+ii) == -9999.99 || *(rap+ii) == -9999.99 ){
      *(jr+ii)= 9999.99;
      continue;
    }
    if ( (*(rap+ii) - *(rperi+ii)) / *(rap+ii) < 0.000001 ){//circular
      *(jr+ii) = 0.;
      continue;
    }
    //Setup function
    (params+tid)->ER= *(ER+ii);
    (params+tid)->Lz22= 0.5 * *(Lz+ii) * *(Lz+ii);
    (JRInt+tid)->function = &JRAdiabaticIntegrand;
    (JRInt+tid)->params = params+tid;
    //Integrate
    *(jr+ii)= gsl_integration_glfixed (JRInt+tid,*(rperi+ii),*(rap+ii),T)
      * sqrt(2.) / M_PI;
  }
  free(JRInt);
  free(params);
  gsl_integration_glfixed_table_free ( T );
}
void calcJzAdiabatic(int ndata,
		     double * jz,
		     double * zmax,
		     double * R,
		     double * Ez,
		     int nargs,
		     struct potentialArg * actionAngleArgs,
		     int order){
  int ii, tid, nthreads;
#ifdef _OPENMP
  nthreads = omp_get_max_threads();
#else
  nthreads = 1;
#endif
  gsl_function * JzInt= (gsl_function *) malloc ( nthreads * sizeof(gsl_function) );
  struct JzAdiabaticArg * params= (struct JzAdiabaticArg *) malloc ( nthreads * sizeof (struct JzAdiabaticArg) );
  for (tid=0; tid < nthreads; tid++){
    (params+tid)->nargs= nargs;
    (params+tid)->actionAngleArgs= actionAngleArgs;
  }
  //Setup integrator
  gsl_integration_glfixed_table * T= gsl_integration_glfixed_table_alloc (order);
  UNUSED int chunk= CHUNKSIZE;
#pragma omp parallel for schedule(static,chunk)				\
  private(tid,ii)							\
  shared(jz,zmax,JzInt,params,T,Ez,R)
  for (ii=0; ii < ndata; ii++){
#ifdef _OPENMP
    tid= omp_get_thread_num();
#else
    tid = 0;
#endif
    if ( *(zmax+ii) == -9999.99 ){
      *(jz+ii)= 9999.99;
      continue;
    }
    if ( *(zmax+ii) < 0.000001 ){//circular
      *(jz+ii) = 0.;
      continue;
    }
    //Setup function
    (params+tid)->Ez= *(Ez+ii);
    (params+tid)->R= *(R+ii);
    (JzInt+tid)->function = &JzAdiabaticIntegrand;
    (JzInt+tid)->params = params+tid;
    //Integrate
    *(jz+ii)= gsl_integration_glfixed (JzInt+tid,0.,*(zmax+ii),T)
      * 2 * sqrt(2.) / M_PI;
  }
  free(JzInt);
  free(params);
  gsl_integration_glfixed_table_free ( T );
}
void calcRapRperi(int ndata,
		  double * rperi,
		  double * rap,
		  double * R,
		  double * ER,
		  double * Lz,
		  int nargs,
		  struct potentialArg * actionAngleArgs){
  int ii, tid, nthreads;
#ifdef _OPENMP
  nthreads = omp_get_max_threads();
#else
  nthreads = 1;
#endif
  double peps, meps;
  gsl_function * JRRoot= (gsl_function *) malloc ( nthreads * sizeof(gsl_function) );
  struct JRAdiabaticArg * params= (struct JRAdiabaticArg *) malloc ( nthreads * sizeof (struct JRAdiabaticArg) );
  //Setup solver
  int status;
  int iter, max_iter = 100;
  const gsl_root_fsolver_type *T;
  double R_lo, R_hi;
  struct pragmasolver *s= (struct pragmasolver *) malloc ( nthreads * sizeof (struct pragmasolver) );
  T = gsl_root_fsolver_brent;
  for (tid=0; tid < nthreads; tid++){
    (params+tid)->nargs= nargs;
    (params+tid)->actionAngleArgs= actionAngleArgs;
    (s+tid)->s= gsl_root_fsolver_alloc (T);
  }
  UNUSED int chunk= CHUNKSIZE;
  gsl_set_error_handler_off();
#pragma omp parallel for schedule(static,chunk)				\
  private(tid,ii,iter,status,R_lo,R_hi,meps,peps)			\
  shared(rperi,rap,JRRoot,params,s,R,ER,Lz,max_iter)
  for (ii=0; ii < ndata; ii++){
#ifdef _OPENMP
    tid= omp_get_thread_num();
#else
    tid = 0;
#endif
    //Setup function
    (params+tid)->ER= *(ER+ii);
    (params+tid)->Lz22= 0.5 * *(Lz+ii) * *(Lz+ii);
    (JRRoot+tid)->params = params+tid;
    (JRRoot+tid)->function = &JRAdiabaticIntegrandSquared;
    //Find starting points for minimum
    if ( fabs(GSL_FN_EVAL(JRRoot+tid,*(R+ii))) < 0.0000001){ //we are at rap or rperi
      peps= GSL_FN_EVAL(JRRoot+tid,*(R+ii)+0.0000001);
      meps= GSL_FN_EVAL(JRRoot+tid,*(R+ii)-0.0000001);
      if ( fabs(peps) < 0.00000001 && fabs(meps) < 0.00000001 && peps*meps >= 0.) {//circular
	*(rperi+ii) = *(R+ii);
	*(rap+ii) = *(R+ii);
      }
      else if ( peps < 0. && meps > 0. ) {//umax
	*(rap+ii)= *(R+ii);
	R_lo= 0.9 * (*(R+ii) - 0.0000001);
	R_hi= *(R+ii) - 0.00000001;
	while ( GSL_FN_EVAL(JRRoot+tid,R_lo) >= 0. && R_lo > 0.000000001){
	  R_hi= R_lo; //this makes sure that brent evaluates using previous
	  R_lo*= 0.9;
	}
	//Find root
	status = gsl_root_fsolver_set ((s+tid)->s, JRRoot+tid, R_lo, R_hi);
	if (status == GSL_EINVAL) {
	  *(rperi+ii) = 0.;//Assume zero if below 0.000000001
	  continue;
	}
	iter= 0;
	do
	  {
	    iter++;
	    status = gsl_root_fsolver_iterate ((s+tid)->s);
	    R_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
	    R_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
	    status = gsl_root_test_interval (R_lo, R_hi,
					     9.9999999999999998e-13,
					     4.4408920985006262e-16);
	  }
	while (status == GSL_CONTINUE && iter < max_iter);
	// LCOV_EXCL_START
	if (status == GSL_EINVAL) {//Shouldn't ever get here
	  *(rperi+ii) = -9999.99;
	  *(rap+ii) = -9999.99;
	  continue;
	}
	// LCOV_EXCL_STOP
	*(rperi+ii) = gsl_root_fsolver_root ((s+tid)->s);
      }
      else if ( peps > 0. && meps < 0. ){//umin
	*(rperi+ii)= *(R+ii);
	R_lo= *(R+ii) + 0.0000001;
	R_hi= 1.1 * (*(R+ii) + 0.0000001);
	while ( GSL_FN_EVAL(JRRoot+tid,R_hi) >= 0. && R_hi < 37.5) {
	  R_lo= R_hi; //this makes sure that brent evaluates using previous
	  R_hi*= 1.1;
	}
	//Find root
	status = gsl_root_fsolver_set ((s+tid)->s, JRRoot+tid, R_lo, R_hi);
	if (status == GSL_EINVAL) {
	  *(rperi+ii) = -9999.99;
	  *(rap+ii) = -9999.99;
	  continue;
	}
	iter= 0;
	do
	  {
	    iter++;
	    status = gsl_root_fsolver_iterate ((s+tid)->s);
	    R_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
	    R_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
	    status = gsl_root_test_interval (R_lo, R_hi,
					     9.9999999999999998e-13,
					     4.4408920985006262e-16);
	  }
	while (status == GSL_CONTINUE && iter < max_iter);
	// LCOV_EXCL_START
	if (status == GSL_EINVAL) {//Shouldn't ever get here
	  *(rperi+ii) = -9999.99;
	  *(rap+ii) = -9999.99;
	  continue;
	}
	// LCOV_EXCL_STOP
	*(rap+ii) = gsl_root_fsolver_root ((s+tid)->s);
      }
    }
    else {
      R_lo= 0.9 * *(R+ii);
      R_hi= *(R+ii);
      while ( GSL_FN_EVAL(JRRoot+tid,R_lo) >= 0. && R_lo > 0.000000001){
	R_hi= R_lo; //this makes sure that brent evaluates using previous
	R_lo*= 0.9;
      }
      R_hi= (R_lo < 0.9 * *(R+ii)) ? R_lo / 0.9 / 0.9: *(R+ii);
      //Find root
      status = gsl_root_fsolver_set ((s+tid)->s, JRRoot+tid, R_lo, R_hi);
      if (status == GSL_EINVAL) {
	*(rperi+ii) = 0.;//Assume zero if below 0.000000001
      } else {
	iter= 0;
	do
	  {
	    iter++;
	    status = gsl_root_fsolver_iterate ((s+tid)->s);
	    R_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
	    R_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
	    status = gsl_root_test_interval (R_lo, R_hi,
					     9.9999999999999998e-13,
					     4.4408920985006262e-16);
	  }
	while (status == GSL_CONTINUE && iter < max_iter);
	// LCOV_EXCL_START
	if (status == GSL_EINVAL) {//Shouldn't ever get here
	  *(rperi+ii) = -9999.99;
	  *(rap+ii) = -9999.99;
	  continue;
	}
	// LCOV_EXCL_STOP
	*(rperi+ii) = gsl_root_fsolver_root ((s+tid)->s);
      }
      //Find starting points for maximum
      R_lo= *(R+ii);
      R_hi= 1.1 * *(R+ii);
      while ( GSL_FN_EVAL(JRRoot+tid,R_hi) > 0. && R_hi < 37.5) {
	R_lo= R_hi; //this makes sure that brent evaluates using previous
	R_hi*= 1.1;
      }
      R_lo= (R_hi > 1.1 * *(R+ii)) ? R_hi / 1.1 / 1.1: *(R+ii);
      //Find root
      status = gsl_root_fsolver_set ((s+tid)->s, JRRoot+tid, R_lo, R_hi);
      if (status == GSL_EINVAL) {
	*(rperi+ii) = -9999.99;
	*(rap+ii) = -9999.99;
	continue;
      }
      iter= 0;
      do
	{
	  iter++;
	  status = gsl_root_fsolver_iterate ((s+tid)->s);
	  R_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
	  R_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
	  status = gsl_root_test_interval (R_lo, R_hi,
					   9.9999999999999998e-13,
					   4.4408920985006262e-16);
	}
      while (status == GSL_CONTINUE && iter < max_iter);
      // LCOV_EXCL_START
      if (status == GSL_EINVAL) {//Shouldn't ever get here
	*(rperi+ii) = -9999.99;
	*(rap+ii) = -9999.99;
	continue;
      }
      // LCOV_EXCL_STOP
      *(rap+ii) = gsl_root_fsolver_root ((s+tid)->s);
    }
  }
  gsl_set_error_handler (NULL);
  for (tid=0; tid < nthreads; tid++)
    gsl_root_fsolver_free( (s+tid)->s);
  free(s);
  free(JRRoot);
  free(params);
}
void calcZmax(int ndata,
	      double * zmax,
	      double * z,
	      double * R,
	      double * Ez,
	      int nargs,
	      struct potentialArg * actionAngleArgs){
  int ii, tid, nthreads;
#ifdef _OPENMP
  nthreads = omp_get_max_threads();
#else
  nthreads = 1;
#endif
  gsl_function * JzRoot= (gsl_function *) malloc ( nthreads * sizeof(gsl_function) );
  struct JzAdiabaticArg * params= (struct JzAdiabaticArg *) malloc ( nthreads * sizeof (struct JzAdiabaticArg) );
  //Setup solver
  int status;
  int iter, max_iter = 100;
  const gsl_root_fsolver_type *T;
  double z_lo, z_hi;
  struct pragmasolver *s= (struct pragmasolver *) malloc ( nthreads * sizeof (struct pragmasolver) );
  T = gsl_root_fsolver_brent;
  for (tid=0; tid < nthreads; tid++){
    (params+tid)->nargs= nargs;
    (params+tid)->actionAngleArgs= actionAngleArgs;
    (s+tid)->s= gsl_root_fsolver_alloc (T);
  }
  UNUSED int chunk= CHUNKSIZE;
  gsl_set_error_handler_off();
#pragma omp parallel for schedule(static,chunk)				\
  private(tid,ii,iter,status,z_lo,z_hi)				\
  shared(zmax,JzRoot,params,s,z,Ez,R,max_iter)
  for (ii=0; ii < ndata; ii++){
#ifdef _OPENMP
    tid= omp_get_thread_num();
#else
    tid = 0;
#endif
    //Setup function
    (params+tid)->Ez= *(Ez+ii);
    (params+tid)->R= *(R+ii);
    (JzRoot+tid)->function = &JzAdiabaticIntegrandSquared;
    (JzRoot+tid)->params = params+tid;
    //Find starting points for minimum
    if ( fabs(GSL_FN_EVAL(JzRoot+tid,*(z+ii))) < 0.0000001){ //we are at zmax
      *(zmax+ii)= fabs( *(z+ii) );
    }
    else {
      z_lo= fabs ( *(z+ii) );
      z_hi= ( *(z+ii) == 0. ) ? 0.1: 1.1 * fabs( *(z+ii) );
      while ( GSL_FN_EVAL(JzRoot+tid,z_hi) >= 0. && z_hi < 37.5) {
	z_lo= z_hi; //this makes sure that brent evaluates using previous
	z_hi*= 1.1;
      }
      //Find root
      status = gsl_root_fsolver_set ((s+tid)->s, JzRoot+tid, z_lo, z_hi);
      if (status == GSL_EINVAL) {
	*(zmax+ii) = -9999.99;
	continue;
      }
      iter= 0;
      do
	{
	  iter++;
	  status = gsl_root_fsolver_iterate ((s+tid)->s);
	  z_lo = gsl_root_fsolver_x_lower ((s+tid)->s);
	  z_hi = gsl_root_fsolver_x_upper ((s+tid)->s);
	  status = gsl_root_test_interval (z_lo, z_hi,
					   9.9999999999999998e-13,
					   4.4408920985006262e-16);
	}
      while (status == GSL_CONTINUE && iter < max_iter);
      // LCOV_EXCL_START
      if (status == GSL_EINVAL) {//Shouldn't ever get here
	*(zmax+ii) = -9999.99;
	continue;
      }
      // LCOV_EXCL_STOP
      *(zmax+ii) = gsl_root_fsolver_root ((s+tid)->s);
    }
  }
  gsl_set_error_handler (NULL);
  for (tid=0; tid < nthreads; tid++)
    gsl_root_fsolver_free( (s+tid)->s);
  free(s);
  free(JzRoot);
  free(params);
}
double JRAdiabaticIntegrand(double R,
			   void * p){
  return sqrt(JRAdiabaticIntegrandSquared(R,p));
}
double JRAdiabaticIntegrandSquared(double R,
				  void * p){
  struct JRAdiabaticArg * params= (struct JRAdiabaticArg *) p;
  return params->ER - evaluatePotentials(R,0.,params->nargs,params->actionAngleArgs) - params->Lz22 / R / R;
}
double JzAdiabaticIntegrand(double z,
			    void * p){
  return sqrt(JzAdiabaticIntegrandSquared(z,p));
}
double JzAdiabaticIntegrandSquared(double z,
				   void * p){
  struct JzAdiabaticArg * params= (struct JzAdiabaticArg *) p;
  return params->Ez - evaluateVerticalPotentials(params->R,z,
						 params->nargs,
						 params->actionAngleArgs);
}
double evaluateVerticalPotentials(double R, double z,
				  int nargs,
				  struct potentialArg * actionAngleArgs){
  return evaluatePotentials(R,z,nargs,actionAngleArgs)
    -evaluatePotentials(R,0.,nargs,actionAngleArgs);
}

max_stars_repo_name: fardal/galpy

max_issues_repo_name: fardal/galpy

max_forks_repo_name: fardal/galpy


text:
//
//  her2k.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 1/6/12.
//  Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//

#ifndef _her2k_h
#define _her2k_h

/// @file her2k.h Performs multiplication of two complex matrices resulting in a Hermitian matrix.

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

namespace LATL
{
   /// @brief Performs multiplcation of complex matrices resulting in a Hermitian matrix.
   ///
   /// For complex matrices A and B, complex Hermitian matrix C, and complex scalar alpha and real
   /// scalar beta
   ///
   ///        C := alpha*A*B'+conj(alpha)*B*A'+beta*C
   /// or
   ///
   ///        C := alpha*A'*B+conj(alpha)*B'*A+beta*C
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param uplo Specifies whether the upper or lower triangular part of the Hermitian matrix C
   /// is to be referenced:
   ///
   ///        if uplo = 'U' or 'u' then C is upper triangular,
   ///        if uplo = 'L' or 'l' then C is lower triangular.
   /// @param trans Specifies the operation to be perfomed as follows:
   ///
   ///        if trans = 'N' or 'n' then C := alpha*A*B'+conj(alpha)*B*A'+beta*C
   ///        if trans = 'C' or 'c' then C := alpha*A'*B+conj(alpha)*B'*A+beta*C
   /// @param n Specifies the order of the complex Hermitian matrix C.  n>=0
   /// @param k Specifies the other dimension of the complex matrices A and B.
   ///
   ///        if trans = 'N' or 'n' then A and B are n-by-k
   ///        if trans = 'C' or 'c' then A and B are k-by-n
   /// @param alpha Complex scalar.
   /// @param A Pointer to complex matrix.
   /// @param ldA Column length of the matrix A. If trans = 'N' or 'n' ldA>=n, otherwise ldA>=k.
   /// @param B Pointer to complex matrix.
   /// @param ldB Column length of the matrix B. If trans = 'N' or 'n' ldB>=n, otherwise ldB>=k.
   /// @param beta Real scalar.
   /// @param C Pointer to complex Hermitian n-by-n matrix C.
   /// Only the upper or lower triangular part of C is referenced, depending on the value of uplo above.
   /// @param ldC Column length of the matrix C.  ldC>=n
   /// @ingroup BLAS

   template <typename real_t>
   int HER2K(char uplo, char trans, int_t n, int_t k, complex<real_t> alpha, complex<real_t> *A, int_t ldA, complex<real_t> *B, int_t ldB, real_t beta, complex<real_t> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      using std::real;
      using std::conj;
      const complex<real_t> zero(0.0,0.0);
      const complex<real_t> one(1.0,0.0);
      int_t i,j,l;
      complex<real_t> *a,*b,*c,*at,*bt;
      complex<real_t> s,t;

      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;
      else if((n==0)||(((alpha==zero)||(k==0))&&(beta==one)))
         return 0;

      if(alpha==zero)
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<j;i++)
                  c[i]*=beta;
               c[j]=beta*real(c[j]);
               c+=ldC;
            }

         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               c[j]=beta*real(c[j]);
               for(i=j+1;i<n;i++)
                  c[i]*=beta;
               c+=ldC;
            }
         }
      }
      else if(trans=='N')
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<j;i++)
                  c[i]*=beta;
               c[j]=beta*real(c[j]);
               a=A;
               b=B;
               for(l=0;l<k;l++)
               {
                  s=alpha*conj(b[j]);
                  t=conj(alpha*a[j]);
                  for(i=0;i<j;i++)
                     c[i]+=s*a[i]+t*b[i];
                  c[j]=real(c[j])+real(a[j]*s+b[j]*t);
                  a+=ldA;
                  b+=ldB;
               }
               c+=ldC;
            }
         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               c[j]=beta*real(c[j]);
               for(i=j+1;i<n;i++)
                  c[i]*=beta;
               a=A;
               b=B;
               for(l=0;l<k;l++)
               {
                  s=alpha*conj(b[j]);
                  t=conj(alpha*a[j]);
                  c[j]=real(c[j])+real(s*a[j]+t*b[j]);
                  for(i=j+1;i<n;i++)
                     c[i]+=s*a[i]+t*b[i];
                  a+=ldA;
                  b+=ldB;
               }
               c+=ldC;
            }
         }
      }
      else
      {
         if(uplo=='U')
         {
            c=C;
            at=A;
            bt=B;
            for(j=0;j<n;j++)
            {
               a=A;
               b=B;
               for(i=0;i<=j;i++)
               {
                  s=zero;
                  t=zero;
                  for(l=0;l<k;l++)
                  {
                     s+=conj(b[l])*at[l];
                     t+=conj(a[l])*bt[l];
                  }
                  if(i<j)
                     c[i]=alpha*t+conj(alpha)*s+beta*c[i];
                  else
                     c[j]=real(alpha*t+conj(alpha)*s)+beta*real(c[j]);
                  a+=ldA;
                  b+=ldB;
               }
               at+=ldA;
               bt+=ldB;
               c+=ldC;
            }
         }
         else
         {
            c=C;
            at=A;
            bt=B;
            for(j=0;j<n;j++)
            {
               a=A+j*ldA;
               b=B+j*ldB;
               for(i=j;i<n;i++)
               {
                  s=zero;
                  t=zero;
                  for(l=0;l<k;l++)
                  {
                     s+=conj(b[l])*at[l];
                     t+=conj(a[l])*bt[l];
                  }
                  if(i>j)
                     c[i]=alpha*t+conj(alpha)*s+beta*c[i];
                  else
                     c[j]=real(alpha*t+conj(alpha)*s)+beta*real(c[j]);
                  a+=ldA;
                  b+=ldB;
               }
               at+=ldA;
               bt+=ldB;
               c+=ldC;
            }
         }
      }
      return 0;
   }

#ifdef __latl_cblas
#include <cblas.h>

   template <> int HER2K<float>(char uplo, char trans, int_t n, int_t k, complex<float> alpha, complex<float> *A, int_t ldA, complex<float> *B, int_t ldB, float beta, complex<float> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_cher2k(CblasColMajor,Uplo,Trans,n,k,&alpha,A,ldA,B,ldB,beta,C,ldC);

      return 0;
   }

   template <> int HER2K<double>(char uplo, char trans, int_t n, int_t k, complex<double> alpha, complex<double> *A, int_t ldA, complex<double> *B, int_t ldB, double beta, complex<double> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldB<((trans=='N')?n:k))
         return -9;
      else if(ldC<n)
         return -12;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_zher2k(CblasColMajor,Uplo,Trans,n,k,&alpha,A,ldA,B,ldB,beta,C,ldC);

      return 0;
   }
#endif

}
#endif

ext: h


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_errno.h>

BC_INLINE2(GSL_ERROR_SELECT_2,int,int,int)
BC_INLINE3(GSL_ERROR_SELECT_3,int,int,int,int)
BC_INLINE4(GSL_ERROR_SELECT_4,int,int,int,int,int)
BC_INLINE5(GSL_ERROR_SELECT_5,int,int,int,int,int,int)
BC_INLINE2VOID(GSL_STATUS_UPDATE,int*,int)


text:
#pragma once

#include "CesiumGltfReader/Library.h"

#include <CesiumAsync/AsyncSystem.h>
#include <CesiumAsync/Future.h>
#include <CesiumAsync/HttpHeaders.h>
#include <CesiumAsync/IAssetAccessor.h>
#include <CesiumGltf/Model.h>
#include <CesiumJsonReader/ExtensionReaderContext.h>
#include <CesiumJsonReader/IExtensionJsonHandler.h>

#include <gsl/span>

#include <functional>
#include <memory>
#include <optional>
#include <string>
#include <vector>

namespace CesiumGltfReader {

/**
 * @brief The result of reading a glTF model with
 * {@link GltfReader::readModel}.
 */
struct CESIUMGLTFREADER_API ModelReaderResult {
  /**
   * @brief The read model, or std::nullopt if the model could not be read.
   */
  std::optional<CesiumGltf::Model> model;

  /**
   * @brief Errors, if any, that occurred during the load process.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warnings, if any, that occurred during the load process.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief The result of reading an image with
 * {@link GltfReader::readImage}.
 */
struct CESIUMGLTFREADER_API ImageReaderResult {

  /**
   * @brief The {@link ImageCesium} that was read.
   *
   * This will be `std::nullopt` if the image could not be read.
   */
  std::optional<CesiumGltf::ImageCesium> image;

  /**
   * @brief Error messages that occurred while trying to read the image.
   */
  std::vector<std::string> errors;

  /**
   * @brief Warning messages that occurred while reading the image.
   */
  std::vector<std::string> warnings;
};

/**
 * @brief Options for how to read a glTF.
 */
struct CESIUMGLTFREADER_API ReadModelOptions {
  /**
   * @brief Whether data URLs in buffers and images should be automatically
   * decoded as part of the load process.
   */
  bool decodeDataUrls = true;

  /**
   * @brief Whether data URLs should be cleared after they are successfully
   * decoded.
   *
   * This reduces the memory usage of the model.
   */
  bool clearDecodedDataUrls = true;

  /**
   * @brief Whether embedded images in {@link Model::buffers} should be
   * automatically decoded as part of the load process.
   *
   * The {@link ImageSpec::mimeType} property is ignored, and instead the
   * [stb_image](https://github.com/nothings/stb) library is used to decode
   * images in `JPG`, `PNG`, `TGA`, `BMP`, `PSD`, `GIF`, `HDR`, or `PIC` format.
   */
  bool decodeEmbeddedImages = true;

  /**
   * @brief Whether geometry compressed using the `KHR_draco_mesh_compression`
   * extension should be automatically decoded as part of the load process.
   */
  bool decodeDraco = true;
};

/**
 * @brief Reads glTF models and images.
 */
class CESIUMGLTFREADER_API GltfReader {
public:
  /**
   * @brief Constructs a new instance.
   */
  GltfReader();

  /**
   * @brief Gets the context used to control how extensions are loaded from glTF
   * files.
   */
  CesiumJsonReader::ExtensionReaderContext& getExtensions();

  /**
   * @brief Gets the context used to control how extensions are loaded from glTF
   * files.
   */
  const CesiumJsonReader::ExtensionReaderContext& getExtensions() const;

  /**
   * @brief Reads a glTF or binary glTF (GLB) from a buffer.
   *
   * @param data The buffer from which to read the glTF.
   * @param options Options for how to read the glTF.
   * @return The result of reading the glTF.
   */
  ModelReaderResult readModel(
      const gsl::span<const std::byte>& data,
      const ReadModelOptions& options = ReadModelOptions()) const;

  /**
   * @brief Accepts the result of {@link readModel} and resolves any remaining
   * external buffers and images.
   *
   * @param asyncSystem The async system to use for resolving external data.
   * @param baseUrl The base url that all the external uris are relative to.
   * @param headers The http headers needed to make any external data requests.
   * @param pAssetAccessor The asset accessor to use to request the external
   * buffers and images.
   * @param result The result of the synchronous readModel invocation.
   */
  static CesiumAsync::Future<ModelReaderResult> resolveExternalData(
      CesiumAsync::AsyncSystem asyncSystem,
      const std::string& baseUrl,
      const CesiumAsync::HttpHeaders& headers,
      std::shared_ptr<CesiumAsync::IAssetAccessor> pAssetAccessor,
      ModelReaderResult&& result);

  /**
   * @brief Reads an image from a buffer.
   *
   * The [stb_image](https://github.com/nothings/stb) library is used to decode
   * images in `JPG`, `PNG`, `TGA`, `BMP`, `PSD`, `GIF`, `HDR`, or `PIC` format.
   *
   * @param data The buffer from which to read the image.
   * @return The result of reading the image.
   */
  static ImageReaderResult readImage(const gsl::span<const std::byte>& data);

private:
  CesiumJsonReader::ExtensionReaderContext _context;
};

} // namespace CesiumGltfReader


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 *  \file smooth.c
 *  \brief smoothing functions using FFT provided by GSL.
 *
 *  References: Press et al., Numerical Recipes, Chapter 13.1.
 *              GSL documents
 */

///////////////THIS FILE IS DEPRECATED///////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <mpi.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>

#include "brains.h"

int nd_fft;

double *resp, *data_fft;
double *resp_cmp, *data_fft_cmp, *data_fft_inverse;

gsl_fft_real_wavetable * real_data, *real_resp;
gsl_fft_complex_wavetable * cmp_data;
gsl_fft_real_workspace * work_data, *work_resp;
gsl_fft_complex_workspace * work_cmp;

/*!
 * This function initiates workspace for FFT.
 */
void smooth_init(int nv, const double *transv)
{
  //nd_fft = parset.n_vel_recon>=n_vel_data? parset.n_vel_recon : n_vel_data;
  nd_fft = nv;
  //printf("%d\n", nd_fft);
  real_data = gsl_fft_real_wavetable_alloc (nd_fft);
  real_resp = gsl_fft_real_wavetable_alloc (nd_fft);
  cmp_data = gsl_fft_complex_wavetable_alloc (nd_fft);

  work_data = gsl_fft_real_workspace_alloc (nd_fft);
  work_resp = gsl_fft_real_workspace_alloc (nd_fft);
  work_cmp = gsl_fft_complex_workspace_alloc (nd_fft);

  resp = malloc(nd_fft * sizeof(double));
  data_fft = malloc(nd_fft * sizeof(double));

  resp_cmp = malloc(nd_fft * 2 * sizeof(double));
  data_fft_cmp = malloc(nd_fft *2 * sizeof(double));
  data_fft_inverse = malloc(nd_fft *2 * sizeof(double));

  // initialize response and its fft.
  int i;
  double sigV, dV, tot;

  sigV = parset.InstRes / VelUnit;
  dV = transv[1] - transv[0];

  /* setup response, whose negective-time part is wrapped around and stored at the right hand*/
  tot = 0.0;
  for (i = 0; i<nd_fft/2; i++)
  {
    resp[i] = 1.0/sqrt(2.0*M_PI)/sigV * exp(-0.5*(i*dV)*(i*dV)/sigV/sigV);
    tot += resp[i];
  }
  for (i = nd_fft-1; i>= nd_fft/2; i--)
  {
    resp[i] = 1.0/sqrt(2.0*M_PI)/sigV * exp(-0.5*((i-nd_fft)*dV)*((i-nd_fft)*dV)/sigV/sigV);
    tot += resp[i];
  }

  /* normalize response */
  for(i=0; i<nd_fft; i++)
  {
    //resp[i] /= (tot * dV);
    resp[i] /= (tot);
  }

  /* FFT of response */
  gsl_fft_real_transform(resp, 1, nd_fft, real_resp, work_resp);
  gsl_fft_halfcomplex_unpack(resp, resp_cmp, 1, nd_fft);
}

/*!
 * This function finalizes FFT.
 */
void smooth_end()
{

  gsl_fft_real_wavetable_free(real_data);
  gsl_fft_real_wavetable_free(real_resp);
  gsl_fft_complex_wavetable_free(cmp_data);

  gsl_fft_real_workspace_free(work_data);
  gsl_fft_real_workspace_free(work_resp);
  gsl_fft_complex_workspace_free(work_cmp);

  free(resp);
  free(data_fft);

  free(resp_cmp);
  free(data_fft_cmp);
  free(data_fft_inverse);

}

/*!
 * This function performs FFT-smoothing to 2d line.
 */
void line_gaussian_smooth_2D_FFT(const double *transv, double *fl2d, int nl, int nv)
{
  int i, j;
  double dV;

  dV = transv[1] - transv[0];

  for(j=0; j<nl; j++)
  {
    memcpy(data_fft, &fl2d[j*nv], nv*sizeof(double));

    /* FFT of line */
    gsl_fft_real_transform(data_fft, 1, nd_fft, real_data, work_data);
    gsl_fft_halfcomplex_unpack(data_fft, data_fft_cmp, 1, nd_fft);

    /* complex multiply and inverse FFT */
    for(i=0; i<nd_fft; i++)
    {
      data_fft_inverse[i*2] = data_fft_cmp[i*2]*resp_cmp[i*2] - data_fft_cmp[i*2+1]*resp_cmp[i*2+1];
      data_fft_inverse[i*2+1] = data_fft_cmp[i*2]*resp_cmp[i*2+1] + data_fft_cmp[i*2+1]*resp_cmp[i*2];
    }
    gsl_fft_complex_inverse(data_fft_inverse, 1, nd_fft, cmp_data, work_cmp);

    for(i=0; i<nv; i++)
    {
      /* take into account the velocity grid width */
      //fl2d[j*nv + i] = data_fft_inverse[i*2] * dV;
      fl2d[j*nv + i] = data_fft_inverse[i*2];
    }
  }
  return;
}

/*!
 * test smoothing function, assumed line profile.
 */
void line_gaussian_smooth_2D_FFT_test(const double *transv, double *fl2d, int nl, int nv)
{
  int i, j;
  double sigV, dV, tot;

  sigV = parset.InstRes / VelUnit;
  dV = transv[1] - transv[0];

  tot = 0.0;
  for (i = 0; i<nd_fft/2; i++)
  {
    resp[i] = 1.0/sqrt(2.0*M_PI)/sigV * exp(-0.5*(i*dV)*(i*dV)/sigV/sigV);
    tot += resp[i];
  }
  for (i = nd_fft-1; i>= nd_fft/2; i--)
  {
    resp[i] = 1.0/sqrt(2.0*M_PI)/sigV * exp(-0.5*((i-nd_fft)*dV)*((i-nd_fft)*dV)/sigV/sigV);
    tot += resp[i];
  }

  for(i=0; i<nd_fft; i++)
  {
    resp[i] /= (tot * dV);
  }
  gsl_fft_real_transform(resp, 1, nd_fft, real_resp, work_resp);
  gsl_fft_halfcomplex_unpack(resp, resp_cmp, 1, nd_fft);

  for(j=0; j<nl; j++)
  {
    //memcpy(data_fft, &fl2d[j*nv], nv*sizeof(double));
    for(i=0; i<nv; i++)
    {
      data_fft[i] =   exp( - 0.5 * pow(transv[i]/(1500.0/VelUnit), 2.0));
    }

    gsl_fft_real_transform(data_fft, 1, nd_fft, real_data, work_data);
    gsl_fft_halfcomplex_unpack(data_fft, data_fft_cmp, 1, nd_fft);

    for(i=0; i<nd_fft; i++)
    {
      data_fft_inverse[i*2] = data_fft_cmp[i*2]*resp_cmp[i*2] - data_fft_cmp[i*2+1]*resp_cmp[i*2+1];
      data_fft_inverse[i*2+1] = data_fft_cmp[i*2]*resp_cmp[i*2+1] + data_fft_cmp[i*2+1]*resp_cmp[i*2];
    }
    gsl_fft_complex_inverse(data_fft_inverse, 1, nd_fft, cmp_data, work_cmp);

    for(i=0; i<nv; i++)
    {
      // take into account the velocity grid width
      fl2d[j*nv + i] = data_fft_inverse[i*2] * dV;
    }
  }
  return;
}

/*!
 * convolution test.
 */
void smooth_test()
{
  int i, n = 205;
  double data[n], resp[n];
  double data_cmp[n*2], resp_cmp[n*2];
  double ans[n*2];

  gsl_fft_real_wavetable * real, *real_resp;
  gsl_fft_complex_wavetable * hc;
  gsl_fft_real_workspace * work, *work_resp;
  gsl_fft_complex_workspace * work_cmp;

  for (i = 0; i < n; i++)
  {
    data[i] = 1.0/sqrt(2.0*M_PI)/6.0 * exp(-0.5*(i-50.0)*(i-50.0)/6.0/6.0);
  }

  for (i = 0; i < n/2; i++)
  {
    resp[i] = 1.0/sqrt(2.0*M_PI)/10.0 * exp(-0.5*(i-0.0)*(i-0.0)/10.0/10.0);
  }
  for (i = n-1; i >= n/2; i--)
  {
    resp[i] = 1.0/sqrt(2.0*M_PI)/10.0 * exp(-0.5*(i-n)*(i-n)/10.0/10.0);
  }

  for (i = 0; i < n; i++)
  {
    printf ("%d: %e\n", i, data[i]);
  }
  printf ("\n");

  work = gsl_fft_real_workspace_alloc (n);
  real = gsl_fft_real_wavetable_alloc (n);

  work_resp = gsl_fft_real_workspace_alloc (n);
  real_resp = gsl_fft_real_wavetable_alloc (n);

  work_cmp = gsl_fft_complex_workspace_alloc(n);
  hc = gsl_fft_complex_wavetable_alloc (n);

  gsl_fft_real_transform (data, 1, n, real, work);
  gsl_fft_real_transform (resp, 1, n, real_resp, work_resp);

  gsl_fft_halfcomplex_unpack(data, data_cmp, 1, n);
  gsl_fft_halfcomplex_unpack(resp, resp_cmp, 1, n);

  for(i=0; i<n; i++)
  {
    ans[i*2] = data_cmp[i*2]*resp_cmp[i*2] - data_cmp[i*2+1]*resp_cmp[i*2+1];
    ans[i*2+1] = data_cmp[i*2]*resp_cmp[i*2+1] + data_cmp[i*2+1]*resp_cmp[i*2];
  }

  gsl_fft_complex_inverse(ans, 1, n, hc, work_cmp);

  for (i = 0; i < n; i++)
  {
    printf ("%d: %e %e\n", i,  ans[i*2], ans[i*2+1]);
  }

  gsl_fft_real_wavetable_free (real);
  gsl_fft_real_wavetable_free (real_resp);
  gsl_fft_real_workspace_free (work);
  gsl_fft_real_workspace_free (work_resp);
  gsl_fft_complex_wavetable_free (hc);
  gsl_fft_complex_workspace_free(work_cmp);
  return;
}

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
#include "del2mat.h"

#define DEL2MAT_MULT   ((void(*)(void))Del2Mat_mult)
#define DEL2MAT_DIAG   ((void(*)(void))Del2Mat_diag)

int main(int argc,char **argv)
{
  PetscInt n;
  PetscScalar h;
  Del2Mat shell;
  Mat A;
  Vec x,b;
  KSP ksp;
  PC  pc;
  /* PETSc initialization  */
  PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
  /* number of nodes in each direction
   * excluding those at the boundary */
  n = 32;
  h = 1.0/(n+1); /* grid spacing */
  /* setup linear system (shell) matrix */
  MatCreate(PETSC_COMM_SELF, &A);
  MatSetSizes(A, n*n*n, n*n*n, n*n*n, n*n*n);
  MatSetType(A, MATSHELL);
  shell.N = n;
  PetscMalloc((n+2)*(n+2)*(n+2)*sizeof(PetscScalar),&shell.F);
  PetscMemzero(shell.F, (n+2)*(n+2)*(n+2)*sizeof(PetscScalar));
  MatShellSetContext(A, (void**)&shell);
  MatShellSetOperation(A, MATOP_MULT,           DEL2MAT_MULT);
  MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, DEL2MAT_MULT);
  MatShellSetOperation(A, MATOP_GET_DIAGONAL,   DEL2MAT_DIAG);
  MatSetUp(A);
  /* setup linear system vectors */
  MatCreateVecs(A, &x, &b);
  VecSet(x, 0);
  VecSet(b, 1);
  /* setup Krylov linear solver */
  KSPCreate(PETSC_COMM_SELF, &ksp);
  KSPGetPC(ksp, &pc);
  KSPSetType(ksp, KSPCG); /* use conjugate gradients */
  PCSetType(pc, PCNONE);  /* with no preconditioning */
  KSPSetFromOptions(ksp);
  /* iteratively solve linear system of equations A*x=b */
  KSPSetOperators(ksp,A,A);
  KSPSolve(ksp, b, x);
  /* scale solution vector to account for grid spacing */
  VecScale(x, h*h);
  /* free memory and destroy objects */
  PetscFree(shell.F);
  VecDestroy(&x);
  VecDestroy(&b);
  MatDestroy(&A);
  KSPDestroy(&ksp);
  /* finalize PETSc */
  PetscFinalize();
  return 0;
}


text:
#ifndef TESTING_COMMON
#define TESTING_COMMON

#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))

// header ------------------------------------------
#include "testing_setting.h"
#undef _IEEE_

#define __STDC_WANT_IEC_60559_TYPES_EXT__
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include <string.h>
#include <stdint.h>
#include <float.h>
#define __STDC_FORMAT_MACROS
#include <inttypes.h>
#include <sys/time.h>
#include <unistd.h>
#include <getopt.h>
#include <ctype.h>
#include <sys/utsname.h>

#if defined (MKL)
//#include <mkl.h>
#include <mkl_cblas.h>
#include <mkl_trans.h>
//#include <mkl_spblas.h> // conflict with Bebop SMC
#else
#include <cblas.h>
#endif

#if defined (ARM)
#define __float128 long double
#endif

#if defined (CUBLAS)
#include <cublas_v2.h>
#elif defined (CUOZBLAS)
#include <cuozblas.h>
#elif defined (OZBLAS)
#include <ozblas.h>
#endif

#if defined (CUDA)
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <cublas_v2.h>
#include <cusparse_v2.h>
#endif

#if defined (CSRMV) || defined (CG)
extern "C" {
#include <bebop/smc/sparse_matrix.h>
#include <bebop/smc/sparse_matrix_ops.h>
#include <bebop/smc/csr_matrix.h>
}
#endif

#if !defined (MPLAPACK)
#define mpreal double
#endif

// ===========================================================
// common for FP_TYPE
#if defined (PREC_Q_D) || defined (PREC_Q_S)
#define PREC_Q
#elif defined (PREC_D_D) || defined (PREC_D_S)
#define PREC_D
#elif defined (PREC_S_S) || defined (PREC_S_D)
#define PREC_S
#elif defined (PREC_32F_PDT_32F)
#define PREC_S
#elif defined (PREC_32F_TC_32F)
#define PREC_S
#elif defined (PREC_32F_TC_32TF)
#define PREC_S
#elif defined (PREC_16F_TC_32F)
#define PREC_H
#elif defined (PREC_16F_PDT_32F)
#define PREC_H
#elif defined (PREC_16F_TC_16F)
#define PREC_H
#elif defined (PREC_16F_PDT_16F)
#define PREC_H
#elif defined (PREC_64F_PDT_64F)
#define PREC_D
#elif defined (PREC_64F_TC_64F)
#define PREC_D
#endif
// ===========================================================

#if defined (MPLAPACK)
#if defined (PREC_Q)
#define MPFR_WANT_FLOAT128
#define _Float128 __float128 // this is needed if MPFR is >= 4.1.0
#endif
#include <mplapack/mpblas_mpfr.h>
#include <mpfr.h> // this is after above
#endif
#undef _Float128

// +++++++++++++++++++++
// FP_TYPE == binary128
// +++++++++++++++++++++
#if defined (PREC_Q)
#define FP_TYPE	__float128
#if defined (ARM)
#define FABS fabs
#define FP_MAX LDBL_MAX
#define FP_MIN LDBL_MIN
#define F_PI 3.1415926535897932384626433832795029L;
#else
#include <quadmath.h>
#define FABS fabsq
#define FP_MAX FLT128_MAX
#define FP_MIN FLT128_MIN
#define F_PI M_PIq
#endif
#if defined (MPLAPACK)
#include <mplapack/mpblas__Float128.h>
#endif

// +++++++++++++++++++++
// FP_TYPE == double-double
// +++++++++++++++++++++
#elif defined (PREC_DD)
#define FP_TYPE dd_real
#define FABS fabs
#if defined (MPLAPACK)
#include <qd/dd_real.h>
#include <mplapack/mpblas_dd.h>
#endif
#define FP_MAX DBL_MAX
#define FP_MIN DBL_MIN
#define F_PI dd_real::_pi

// +++++++++++++++++++++
// FP_TYPE == binary64
// +++++++++++++++++++++
#elif defined (PREC_D)
#define FP_TYPE	double
#define FABS fabs
#define FP_MAX DBL_MAX
#define FP_MIN DBL_MIN
#define F_PI M_PI

// +++++++++++++++++++++
// FP_TYPE == binary32
// +++++++++++++++++++++
#elif defined (PREC_S)
#define FP_TYPE	float
#define FABS fabs
#define FP_MAX FLT_MAX
#define FP_MIN FLT_MIN
#define F_PI M_PI

// +++++++++++++++++++++
// FP_TYPE == binary16
// +++++++++++++++++++++
#elif defined (PREC_H)
#include "half-2.1.0/half.hpp"
using half_float::half;
#define FP_TYPE	half
#define FABS fabs
#define FP_MAX 65504.
#define FP_MIN (2.e-14)-1
#define F_PI M_PI

#endif

// ===========================================================
// reference BLAS
#if defined (MPLAPACK)
#define refRdot												Rdot
#define refRgemv(tran,m,n,al,A,lda,X,ix,bt,Y,iy)			Rgemv(&tran,m,n,al,A,lda,X,ix,bt,Y,iy)
#define refRgemm(tranA,tranB,m,n,k,al,A,lda,B,ldb,bt,C,ldc) Rgemm(&tranA,&tranB,m,n,k,al,A,lda,B,ldb,bt,C,ldc)
#else
#define refRdot	cblas_ddot
#define refRgemv(tran,m,n,al,A,lda,X,ix,bt,Y,iy)			cblas_dgemv(CblasColMajor,ToCblasOp(tran),m,n,al,A,lda,X,ix,bt,Y,iy)
#define refRgemm(tranA,tranB,m,n,k,al,A,lda,B,ldb,bt,C,ldc) cblas_dgemm(CblasColMajor,ToCblasOp(tranA),ToCblasOp(tranB),m,n,k,al,A,lda,B,ldb,bt,C,ldc)
#endif

// ===========================================================
// argment transformation
#if defined (OZBLAS)
#define trgRdot(ha,n,x,ix,y,iy)								rdot(&ha,n,x,ix,y,iy)
#define trgRnrm2(ha,n,x,ix)									rnrm2(&ha,n,x,ix)
#define trgRaxpy(ha,n,a,x,ix,y,iy)							raxpy(&ha,n,a,x,ix,y,iy)
#define	trgRgemv(ha,ta,m,n,al,a,la,x,ix,bt,y,iy)			rgemv(&ha,ta,m,n,al,a,la,x,ix,bt,y,iy)
#define trgRgemm(ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)	rgemm(&ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)
#define	trgRcsrmv(ha,ta,m,n,nnz,al,da,a,ci,rp,x,bt,y)		rcsrmv(&ha,ta,m,n,nnz,al,da,a,ci,rp,x,bt,y)
#define	trgRcg(ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)			rcg(&ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)
#elif defined (CUOZBLAS)
#define trgRdot(ha,n,x,ix,y,iy,r)							rdot(&ha,n,x,ix,y,iy,r)
#define trgRnrm2(ha,n,x,ix,r)								rnrm2(&ha,n,x,ix,r)
#define trgRaxpy(ha,n,a,x,ix,y,iy)							raxpy(&ha,n,a,x,ix,y,iy)
#define	trgRgemv(ha,ta,m,n,al,a,la,x,ix,bt,y,iy)			rgemv(&ha,ta,m,n,al,a,la,x,ix,bt,y,iy)
#define trgRgemm(ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)	rgemm(&ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)
#define	trgRcsrmv(ha,ta,m,n,nnz,al,da,a,ci,rp,x,bt,y)		rcsrmv(&ha,ta,m,n,nnz,al,da,a,ci,rp,x,bt,y)
#define	trgRcg(ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)			rcg(&ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)
#elif defined (CUBLAS) && defined (GEMMEX)
#define trgRgemm(ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)	cublasGemmEx(ha,ToCublasOp(ta),ToCublasOp(tb),m,n,k,&al,A,DATA_TYPE_A,lda,B,DATA_TYPE_B,ldb,&bt,C,DATA_TYPE_C,ldc,COMP_TYPE,ALGO_TYPE)
#elif defined (CUBLAS)
#define trgRdot(ha,n,x,ix,y,iy,r)							rdot(ha,n,x,ix,y,iy,r)
#define trgRnrm2(ha,n,x,ix)									rnrm2(ha,n,x,ix)
#define trgRaxpy(ha,n,a,x,ix,y,iy)							raxpy(ha,n,a,x,ix,y,iy)
#define	trgRgemv(ha,ta,m,n,al,a,la,x,ix,bt,y,iy)			rgemv(ha,ToCublasOp(ta),m,n,&al,a,la,x,ix,&bt,y,iy)
#define trgRgemm(ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)	rgemm(ha,ToCublasOp(ta),ToCublasOp(tb),m,n,k,&al,A,lda,B,ldb,&bt,C,ldc)
#define	trgRcsrmv(ha,ta,m,n,nnz,al,da,a,ci,rp,x,bt,y)		rcsrmv(ha,ToCusparseOp(ta),m,n,nnz,al,da,a,ci,rp,x,bt,y)
//#define	trgRcg(ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)			rcg(&ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)
#elif defined (MPBLAS)
#define trgRdot(ha,n,x,ix,y,iy)								rdot(n,x,ix,y,iy)
#define trgRnrm2(ha,n,x,ix)									rnrm2(n,x,ix)
#define trgRaxpy(ha,n,a,x,ix,y,iy)							raxpy(ha,n,a,x,ix,y,iy)
#define	trgRgemv(ha,ta,m,n,al,a,la,x,ix,bt,y,iy)			rgemv(&ta,m,n,al,a,la,x,ix,bt,y,iy)
#define trgRgemm(ha,ta,tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)	rgemm(&ta,&tb,m,n,k,al,A,lda,B,ldb,bt,C,ldc)
//#define	trgRcg(ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)			rcg(&ha,ta,n,nnz,da,a,ci,rp,b,x,itr,tol)
#endif

// ===========================================================
#if defined (CUOZBLAS)
#define ozblasHandle_t		cuozblasHandle_t
#define ozblasCreate 		cuozblasCreate
#define ozblasDestroy		cuozblasDestroy
#endif

// ================================================
#if defined (PREC_Q_D)
#if defined (CUOZBLAS)
#define rcg				cuozblasRcg<__float128,double>
#define rdot			cuozblasRdot<__float128,double>
#define rnrm2			cuozblasRnrm2<__float128,double>
#define raxpy			cuozblasRaxpy<__float128>
#define rgemv			cuozblasRgemv<__float128,double>
#define rgemm			cuozblasRgemm<__float128,double>
#define rcsrmv			cuozblasRcsrmv<__float128,double>
#define rcsrmvSplitA	cuozblasRcsrmvSplitA<__float128,double>
#elif defined (OZBLAS)
#define rcg				ozblasRcg<__float128,double>
#define rdot			ozblasRdot<__float128,double>
#define rnrm2			ozblasRnrm2<__float128,double>
#define rgemv			ozblasRgemv<__float128,double>
#define rgemm			ozblasRgemm<__float128,double>
#define rcsrmv			ozblasRcsrmv<__float128,double>
#define rcsrmvSplitA	ozblasRcsrmvSplitA<__float128,double>
#elif defined (MPBLAS)
#define	rdot			Rdot
#define	rgemv			Rgemv
#define	rgemm			Rgemm
#endif

// ================================================
#elif defined (PREC_Q_S)
#if defined (CUOZBLAS)
#define rcg				cuozblasRcg<__float128,float>
#define rdot			cuozblasRdot<__float128,float>
#define rnrm2			cuozblasRnrm2<__float128,float>
#define raxpy			cuozblasRaxpy<__float128>
#define rgemv			cuozblasRgemv<__float128,float>
#define rgemm			cuozblasRgemm<__float128,float>
#define rcsrmv			cuozblasRcsrmv<__float128,float>
#define rcsrmvSplitA	cuozblasRcsrmvSplitA<__float128,float>
#elif defined (OZBLAS)
#define rcg				ozblasRcg<__float128,float>
#define rdot			ozblasRdot<__float128,float>
#define rnrm2			ozblasRnrm2<__float128,float>
#define rgemv			ozblasRgemv<__float128,float>
#define rgemm			ozblasRgemm<__float128,float>
#define rcsrmv			ozblasRcsrmv<__float128,float>
#define rcsrmvSplitA	ozblasRcsrmvSplitA<__float128,float>
#endif

// ================================================
#elif defined (PREC_D_S)
#if defined (CUOZBLAS)
#define rcg				cuozblasRcg<double,float>
#define rdot			cuozblasRdot<double,float>
#define rnrm2			cuozblasRnrm2<double,float>
#define raxpy			cuozblasRaxpy<double>
#define rgemv			cuozblasRgemv<double,float>
#define rgemm			cuozblasRgemm<double,float>
#define rcsrmv			cuozblasRcsrmv<double,float>
#define rcsrmvSplitA	cuozblasRcsrmvSplitA<double,float>
#elif defined (OZBLAS)
#define rcg				ozblasRcg<double,float>
#define rdot			ozblasRdot<double,float>
#define rnrm2			ozblasRnrm2<double,float>
#define rgemv			ozblasRgemv<double,float>
#define rgemm			ozblasRgemm<double,float>
#define rcsrmv			ozblasRcsrmv<double,float>
#define rcsrmvSplitA	ozblasRcsrmvSplitA<double,float>
#endif

// ================================================
#elif defined (PREC_D_D)
#if defined (CUOZBLAS)
#define rcg				cuozblasRcg<double,double>
#define rdot			cuozblasRdot<double,double>
#define rnrm2			cuozblasRnrm2<double,double>
#define raxpy			cuozblasRaxpy<double>
#define rgemv			cuozblasRgemv<double,double>
//#define rgemm			cuozblasDgemm
#define rgemm			cuozblasRgemm<double,double>
#define rcsrmv			cuozblasRcsrmv<double,double>
#define rcsrmvSplitA	cuozblasRcsrmvSplitA<double,double>
#elif defined (OZBLAS)
#define rcg				ozblasRcg<double,double>
#define rdot			ozblasRdot<double,double>
#define rnrm2			ozblasRnrm2<double,double>
#define rgemv			ozblasRgemv<double,double>
//#define rgemv			ozblasDgemv
#define rgemm			ozblasRgemm<double,double>
#define rcsrmv			ozblasRcsrmv<double,double>
#define rcsrmvSplitA	ozblasRcsrmvSplitA<double,double>
#elif defined (CUBLAS)
//#define rcg			cublasDcg
#define rdot			cublasDdot
#define rnrm2			cublasDnrm2
#define rgemv			cublasDgemv
#define rgemm			cublasDgemm
//#define rcsrmv		cublasDcsrmv
#endif

// ================================================
#elif defined (PREC_S_S)
#if defined (CUOZBLAS)
#define rcg				cuozblasRcg<float,float>
#define rdot			cuozblasRdot<float,float>
#define rnrm2			cuozblasRnrm2<float,float>
#define raxpy			cuozblasRaxpy<float>
#define rgemv			cuozblasRgemv<float,float>
#define rgemm			cuozblasRgemm<float,float>
#define rcsrmv			cuozblasRcsrmv<float,float>
#define rcsrmvSplitA	cuozblasRcsrmvSplitA<float,float>
#elif defined (OZBLAS)
#define rcg				ozblasRcg<float,float>
#define rdot			ozblasRdot<float,float>
#define rnrm2			ozblasRnrm2<float,float>
#define rgemv			ozblasRgemv<float,float>
#define rgemm			ozblasRgemm<float,float>
#define rcsrmv			ozblasRcsrmv<float,float>
#define rcsrmvSplitA	ozblasRcsrmvSplitA<float,float>
#elif defined (CUBLAS)
//#define rcg			cublasScg
#define rdot			cublasSdot
#define rnrm2			cublasSnrm2
#define rgemv			cublasSgemv
#define rgemm			cublasSgemm
//#define rcsrmv		cublasScsrmv
#endif

// ================================================
#elif defined (PREC_S_D)
#if defined (CUOZBLAS)
#define rcg				cuozblasRcg<float,double>
#define rdot			cuozblasRdot<float,double>
#define rnrm2			cuozblasRnrm2<float,double>
#define raxpy			cuozblasRaxpy<float>
#define rgemv			cuozblasRgemv<float,double>
#define rgemm			cuozblasRgemm<float,double>
#define rcsrmv			cuozblasRcsrmv<float,double>
#define rcsrmvSplitA	cuozblasRcsrmvSplitA<float,double>
#elif defined (OZBLAS)
#define rcg				ozblasRcg<float,double>
#define rdot			ozblasRdot<float,double>
#define rnrm2			ozblasRnrm2<float,double>
#define rgemv			ozblasRgemv<float,double>
#define rgemm			ozblasRgemm<float,double>
#define rcsrmv			ozblasRcsrmv<float,double>
#define rcsrmvSplitA	ozblasRcsrmvSplitA<float,double>
#endif

// ================================================
#elif defined (PREC_DD)
#if defined (MPBLAS)
#define	rdot			Rdot
#define	rgemv			Rgemv
#define	rgemm			Rgemm
#endif

#endif

// ================================================
// GemmEx
#if defined (CUBLAS) && defined (GEMMEX)
#define rgemm			cublasGemmEx
// ----------------------------
#if defined (PREC_32F_PDT_32F)
#define DATA_TYPE_A CUDA_R_32F
#define DATA_TYPE_B CUDA_R_32F
#define DATA_TYPE_C CUDA_R_32F
#define COMP_TYPE CUBLAS_COMPUTE_32F_PEDANTIC
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_32F_TC_32F)
#define DATA_TYPE_A CUDA_R_32F
#define DATA_TYPE_B CUDA_R_32F
#define DATA_TYPE_C CUDA_R_32F
#define COMP_TYPE CUBLAS_COMPUTE_32F
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_32F_TC_32TF)
#define DATA_TYPE_A CUDA_R_32F
#define DATA_TYPE_B CUDA_R_32F
#define DATA_TYPE_C CUDA_R_32F
#define COMP_TYPE CUBLAS_COMPUTE_32F_FAST_TF32
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_16F_TC_32F)
#define DATA_TYPE_A CUDA_R_16F
#define DATA_TYPE_B CUDA_R_16F
#define DATA_TYPE_C CUDA_R_16F
#define COMP_TYPE CUBLAS_COMPUTE_32F
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_16F_PDT_32F)
#define DATA_TYPE_A CUDA_R_16F
#define DATA_TYPE_B CUDA_R_16F
#define DATA_TYPE_C CUDA_R_16F
#define COMP_TYPE CUBLAS_COMPUTE_32F_PEDANTIC
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_16F_TC_16F)
#define DATA_TYPE_A CUDA_R_16F
#define DATA_TYPE_B CUDA_R_16F
#define DATA_TYPE_C CUDA_R_16F
#define COMP_TYPE CUBLAS_COMPUTE_16F
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_16F_PDT_16F)
#define DATA_TYPE_A CUDA_R_16F
#define DATA_TYPE_B CUDA_R_16F
#define DATA_TYPE_C CUDA_R_16F
#define COMP_TYPE CUBLAS_COMPUTE_16F_PEDANTIC
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_64F_PDT_64F)
#define DATA_TYPE_A CUDA_R_64F
#define DATA_TYPE_B CUDA_R_64F
#define DATA_TYPE_C CUDA_R_64F
#define COMP_TYPE CUBLAS_COMPUTE_64F_PEDANTIC
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#elif defined (PREC_64F_TC_64F)
#define DATA_TYPE_A CUDA_R_64F
#define DATA_TYPE_B CUDA_R_64F
#define DATA_TYPE_C CUDA_R_64F
#define COMP_TYPE CUBLAS_COMPUTE_64F
#define ALGO_TYPE CUBLAS_GEMM_DEFAULT
// ----------------------------
#endif
#endif

#endif

max_stars_repo_name: wsmoses/Riken-blas

max_issues_repo_name: wsmoses/Riken-blas

max_forks_repo_name: wsmoses/Riken-blas


text:
/**
 *
 * @file testing_sgemm.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Emmanuel Agullo
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:19 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_smain.h"

#undef COMPLEX
#define REAL

int testing_slange(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 3) {
        USAGE("LANGE", "M N LDA",
              "   - M      : number of rows of matrices A and C\n"
              "   - N      : number of columns of matrices B and C\n"
              "   - LDA    : leading dimension of matrix A\n");
        return -1;
    }

    int M     = atoi(argv[0]);
    int N     = atoi(argv[1]);
    int LDA   = atoi(argv[2]);
    int LDAxN = LDA*N;
    int n, u;
    float eps;

    float *A    = (float *)malloc(LDAxN*sizeof(float));
    float             *work = (float*) malloc(max(M,N)*sizeof(float));
    float normplasma, normlapack, result;

    eps = LAPACKE_slamch_work('e');

    printf("\n");
    printf("------ TESTS FOR PLASMA SLANGE ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", M, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 10.\n");

    /*----------------------------------------------------------
     *  TESTING SLANGE
     */

    /* Initialize A, B, C */
    PLASMA_splrnt( M, N, A, LDA, 3436 );

    /* PLASMA SLANGE */
    for(n=0; n<4; n++) {
        normplasma = PLASMA_slange(norm[n], M, N, A, LDA);
        normlapack = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(norm[n]), M, N, A, LDA, work);

        printf("Lapack %e, Plasma %e\n", normlapack, normplasma);

        result = fabs(normplasma - normlapack) / (normlapack * eps);
        switch(norm[n]) {
        case PlasmaMaxNorm:
            /* result should be perfectly equal */
            break;
        case PlasmaInfNorm:
            /* Sum order on the line can differ */
            result = result / (float)N;
            break;
        case PlasmaOneNorm:
            /* Sum order on the column can differ */
            result = result / (float)M;
            break;
        case PlasmaFrobeniusNorm:
            /* Sum oreder on every element can differ */
            result = result / ((float)M * (float)N);
            break;
        }

        printf("***************************************************\n");
        if ( result < 1. ) {
            printf(" ---- TESTING SLANGE (%s)............... PASSED !\n", normstr[n]);
        }
        else {
            printf(" - TESTING SLANGE (%s)... FAILED !\n", normstr[n]);
        }
        printf("***************************************************\n");
    }

    /* Don't perform real tests while lapacke is not correct */
#ifdef COMPLEX
    /* PLASMA SLANTR */
    for(n=0; n<4; n++) {
        for(u=0; u<2; u++) {
            int d;
            for(d=0; d<2; d++) {
                normplasma = PLASMA_slantr(norm[n], uplo[u], diag[d], M, N, A, LDA);
                normlapack = LAPACKE_slantr_work(LAPACK_COL_MAJOR, lapack_const(norm[n]), lapack_const(uplo[u]),
                                                 lapack_const(diag[d]), M, N, A, LDA, work);

                printf("Lapack %e, Plasma %e\n", normlapack, normplasma);

                result = fabs(normplasma - normlapack) / (normlapack * eps);
                switch(norm[n]) {
                case PlasmaMaxNorm:
                    /* result should be perfectly equal */
                    break;
                case PlasmaInfNorm:
                    /* Sum order on the line can differ */
                    result = result / (float)N;
                    break;
                case PlasmaOneNorm:
                    /* Sum order on the column can differ */
                    result = result / (float)M;
                    break;
                case PlasmaFrobeniusNorm:
                    /* Sum oreder on every element can differ */
                    result = result / ((float)M * (float)N);
                    break;
                }

                printf("***************************************************\n");
                if ( result < 1. ) {
                    printf(" ---- TESTING SLANTR (%s, %s, %s)......... PASSED !\n",
                           normstr[n], uplostr[u], diagstr[d]);
                }
                else {
                    printf(" - TESTING SLANTR (%s, %s, %s)... FAILED !\n",
                           normstr[n], uplostr[u], diagstr[d]);
                }
                printf("***************************************************\n");
            }
        }
    }
#endif

    /* PLASMA SLANSY */
    for(n=0; n<4; n++) {
        for(u=0; u<2; u++) {
            normplasma = PLASMA_slansy(norm[n], uplo[u], min(M,N), A, LDA);
            normlapack = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(norm[n]), lapack_const(uplo[u]), min(M,N), A, LDA, work);

            printf("Lapack %e, Plasma %e\n", normlapack, normplasma);

            result = fabs(normplasma - normlapack) / (normlapack * eps);
            switch(norm[n]) {
            case PlasmaMaxNorm:
                /* result should be perfectly equal */
                break;
            case PlasmaInfNorm:
                /* Sum order on the line can differ */
                result = result / (float)N;
                break;
            case PlasmaOneNorm:
                /* Sum order on the column can differ */
                result = result / (float)M;
                break;
            case PlasmaFrobeniusNorm:
                /* Sum oreder on every element can differ */
                result = result / ((float)M * (float)N);
                break;
            }

            printf("***************************************************\n");
            if ( result < 1. ) {
                printf(" ---- TESTING SLANSY (%s, %s)......... PASSED !\n", normstr[n], uplostr[u]);
            }
            else {
                printf(" - TESTING SLANSY (%s, %s)... FAILED !\n", normstr[n], uplostr[u]);
            }
            printf("***************************************************\n");
        }
    }

#ifdef COMPLEX
    /* PLASMA SLANSY */
    {
      int j;
      for (j=0; j<min(M,N); j++) {
        A[j*LDA+j] -= I*cimagf(A[j*LDA+j]);
      }
    }

    for(n=0; n<4; n++) {
        for(u=0; u<2; u++) {
            normplasma = PLASMA_slansy(norm[n], uplo[u], min(M,N), A, LDA);
            normlapack = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(norm[n]), lapack_const(uplo[u]), min(M,N), A, LDA, work);

            printf("Lapack %e, Plasma %e\n", normlapack, normplasma);

            result = fabs(normplasma - normlapack) / (normlapack * eps);
            switch(norm[n]) {
            case PlasmaMaxNorm:
                /* result should be perfectly equal */
                break;
            case PlasmaInfNorm:
                /* Sum order on the line can differ */
                result = result / (float)N;
                break;
            case PlasmaOneNorm:
                /* Sum order on the column can differ */
                result = result / (float)M;
                break;
            case PlasmaFrobeniusNorm:
                /* Sum oreder on every element can differ */
                result = result / ((float)M * (float)N);
                break;
            }

            printf("***************************************************\n");
            if ( result < 1. ) {
                printf(" ---- TESTING SLANSY (%s, %s)......... PASSED !\n", normstr[n], uplostr[u]);
            }
            else {
                printf(" - TESTING SLANSY (%s, %s)... FAILED !\n", normstr[n], uplostr[u]);
            }
            printf("***************************************************\n");
        }
    }
#endif

    free(A); free(work);
    return 0;
}