sc/source/core/opencl/op_addin.cxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  */
   9
  10 #include "op_addin.hxx"
  11
  12 #include <sstream>
  13
  14 namespace sc::opencl {
  15
  16 void OpBesselj::GenSlidingWindowFunction(outputstream &ss,
  17     const std::string &sSymName, SubArguments &vSubArguments)
  18 {
  19     CHECK_PARAMETER_COUNT( 2, 2 );
  20     GenerateFunctionDeclaration( sSymName, vSubArguments, ss );
  21     ss << "{\n";
  22     ss << "    int gid0 = get_global_id(0);\n";
  23     GenerateArg( "x", 0, vSubArguments, ss );
  24     GenerateArg( "N", 1, vSubArguments, ss );
  25     ss << "    double f_2_DIV_PI = 2.0 / M_PI;\n";
  26     ss << "    if( N < 0.0 )\n";
  27     ss << "        return CreateDoubleError(IllegalArgument);\n";
  28     ss << "    if (x == 0.0)\n";
  29     ss << "        return (N == 0.0) ? 1.0 : 0.0;\n";
  30     ss << "    double fSign = ((int)N % 2 == 1 && x < 0.0) ? -1.0 : 1.0;\n";
  31     ss << "    double fX = fabs(x);\n";
  32     ss << "    double fMaxIteration = 9000000.0;\n";
  33     ss << "    double fEstimateIteration = fX * 1.5 + N;\n";
  34     ss << "    bool bAsymptoticPossible = pow(fX,0.4) > N;\n";
  35     ss << "    if (fEstimateIteration > fMaxIteration)\n";
  36     ss << "    {\n";
  37     ss << "        if (bAsymptoticPossible)\n";
  38     ss << "            return fSign * sqrt(f_2_DIV_PI/fX)";
  39     ss << "* cos(fX-N*M_PI_2-M_PI_4);\n";
  40     ss << "        else\n";
  41     ss << "            return CreateDoubleError(NoConvergence);\n";
  42     ss << "    }\n";
  43     ss << "    double epsilon = 1.0e-15;\n";
  44     ss << "    bool bHasfound = false;\n";
  45     ss << "    double k= 0.0;\n";
  46     ss << "    double  u ;\n";
  47     ss << "    double m_bar;\n";
  48     ss << "    double g_bar;\n";
  49     ss << "    double g_bar_delta_u;\n";
  50     ss << "    double g = 0.0;\n";
  51     ss << "    double delta_u = 0.0;\n";
  52     ss << "    double f_bar = -1.0;\n";
  53     ss << "    if (N==0)\n";
  54     ss << "    {\n";
  55     ss << "        u = 1.0;\n";
  56     ss << "        g_bar_delta_u = 0.0;\n";
  57     ss << "        g_bar = - 2.0/fX; \n";
  58     ss << "        delta_u = g_bar_delta_u / g_bar;\n";
  59     ss << "        u = u + delta_u ;\n";
  60     ss << "        g = -1.0 / g_bar; \n";
  61     ss << "        f_bar = f_bar * g;\n";
  62     ss << "        k = 2.0;\n";
  63     ss << "    }\n";
  64     ss << "    if (N!=0)\n";
  65     ss << "    {\n";
  66     ss << "        u=0.0;\n";
  67     ss << "        for (k =1.0; k<= N-1; k = k + 1.0)\n";
  68     ss << "        {\n";
  69     ss << "            m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
  70     ss << "            g_bar_delta_u = - g * delta_u - m_bar * u;\n";
  71     ss << "            g_bar = m_bar - 2.0*k/fX + g;\n";
  72     ss << "            delta_u = g_bar_delta_u / g_bar;\n";
  73     ss << "            u = u + delta_u;\n";
  74     ss << "            g = -1.0/g_bar;\n";
  75     ss << "            f_bar=f_bar * g;\n";
  76     ss << "        }\n";
  77     ss << "        m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
  78     ss << "        g_bar_delta_u = f_bar - g * delta_u - m_bar * u;\n";
  79     ss << "        g_bar = m_bar - 2.0*k/fX + g;\n";
  80     ss << "        delta_u = g_bar_delta_u / g_bar;\n";
  81     ss << "        u = u + delta_u;\n";
  82     ss << "        g = -1.0/g_bar;\n";
  83     ss << "        f_bar = f_bar * g;\n";
  84     ss << "        k = k + 1.0;\n";
  85     ss << "    }\n";
  86     ss << "    do\n";
  87     ss << "    {\n";
  88     ss << "        m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar;\n";
  89     ss << "        g_bar_delta_u = - g * delta_u - m_bar * u;\n";
  90     ss << "        g_bar = m_bar - 2.0*k/fX + g;\n";
  91     ss << "        delta_u = g_bar_delta_u / g_bar;\n";
  92     ss << "        u = u + delta_u;\n";
  93     ss << "        g = -pow(g_bar,-1.0);\n";
  94     ss << "        f_bar = f_bar * g;\n";
  95     ss << "        bHasfound = (fabs(delta_u)<=fabs(u)*epsilon);\n";
  96     ss << "        k = k + 1.0;\n";
  97     ss << "    }\n";
  98     ss << "    while (!bHasfound && k <= fMaxIteration);\n";
  99     ss << "    if (bHasfound)\n";
 100     ss << "        return u * fSign;\n";
 101     ss << "    else\n";
 102     ss << "        return CreateDoubleError(NoConvergence);\n";
 103     ss << "}";
 104 }
 105 void OpGestep::GenSlidingWindowFunction(
 106     outputstream &ss,const std::string &sSymName,
 107     SubArguments &vSubArguments)
 108 {
 109     CHECK_PARAMETER_COUNT( 2, 2 );
 110     GenerateFunctionDeclaration( sSymName, vSubArguments, ss );
 111     ss << "{\n";
 112     ss << "    double tmp=0;\n";
 113     ss << "    int gid0=get_global_id(0);\n";
 114     ss <<"\n";
 115     GenerateArg( "tmp0", 0, vSubArguments, ss );
 116     GenerateArg( "tmp1", 1, vSubArguments, ss );
 117     ss << "    tmp =tmp0 >= tmp1 ? 1 : 0;\n";
 118     ss << "    return tmp;\n";
 119     ss << "}\n";
 120 }
 121
 122 }
 123
 124 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */