LCOV - code coverage report
Current view: top level - src/grid/dgemm - grid_dgemm_prepare_pab.c (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 41 754 5.4 %
Date: 2024-12-12 07:20:04 Functions: 3 12 25.0 %

          Line data    Source code
       1             : /*----------------------------------------------------------------------------*/
       2             : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3             : /*  Copyright 2000-2024 CP2K developers group <https://cp2k.org>              */
       4             : /*                                                                            */
       5             : /*  SPDX-License-Identifier: BSD-3-Clause                                     */
       6             : /*----------------------------------------------------------------------------*/
       7             : 
       8             : #include "grid_dgemm_prepare_pab.h"
       9             : 
      10             : #include <assert.h>
      11             : #include <stdbool.h>
      12             : 
      13             : #include "../common/grid_common.h"
      14             : #include "../common/grid_constants.h"
      15             : #include "grid_dgemm_utils.h"
      16             : 
      17             : struct pab_computation_struct_ {
      18             :   int offset[2];
      19             :   int lmax[2];
      20             :   int lmin[2];
      21             :   double zeta[2];
      22             :   tensor *pab;
      23             :   tensor *pab_prep;
      24             :   int dir1, dir2;
      25             : };
      26             : 
      27             : // *****************************************************************************
      28       44389 : static void grid_prepare_pab_AB(struct pab_computation_struct_ *const tp) {
      29      117327 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
      30      195459 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
      31      300574 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
      32      437913 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
      33      259860 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
      34      562143 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
      35      302283 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
      36      655153 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
      37      352870 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
      38      352870 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
      39      352870 :               idx2(tp->pab_prep[0], coset(lxb, lyb, lzb),
      40      352870 :                    coset(lxa, lya, lza)) = idx2(tp->pab[0], jco, ico);
      41             :             }
      42             :           }
      43             :         }
      44             :       }
      45             :     }
      46             :   }
      47       44389 : }
      48             : 
      49             : // *****************************************************************************
      50           0 : static void grid_prepare_pab_DADB(struct pab_computation_struct_ *const tp) {
      51             :   // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
      52             :   // is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b)
      53             :   // (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx
      54             :   // pgf_{b-1x} - 2*tp->zeta[1]*pgf_{b+1x})
      55             : 
      56           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
      57           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
      58           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
      59           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
      60           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
      61           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
      62           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
      63           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
      64           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
      65           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
      66           0 :               const double pab = idx2(tp->pab[0], jco, ico);
      67           0 :               int ico_l, jco_l;
      68             :               // x  (all safe if lxa = 0, as the spurious added terms have zero
      69             :               // prefactor)
      70             : 
      71           0 :               ico_l = coset(imax(lxa - 1, 0), lya, lza);
      72           0 :               jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
      73           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) += 0.5 * lxa * lxb * pab;
      74             : 
      75           0 :               ico_l = coset(imax(lxa - 1, 0), lya, lza);
      76           0 :               jco_l = coset((lxb + 1), lyb, lzb);
      77           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * tp->zeta[1] * pab;
      78             : 
      79           0 :               ico_l = coset((lxa + 1), lya, lza);
      80           0 :               jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
      81           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) -= tp->zeta[0] * lxb * pab;
      82             : 
      83           0 :               ico_l = coset((lxa + 1), lya, lza);
      84           0 :               jco_l = coset((lxb + 1), lyb, lzb);
      85           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) +=
      86           0 :                   2.0 * tp->zeta[0] * tp->zeta[1] * pab;
      87             : 
      88             :               // y
      89             : 
      90           0 :               ico_l = coset(lxa, imax(lya - 1, 0), lza);
      91           0 :               jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
      92           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) += 0.5 * lya * lyb * pab;
      93             : 
      94           0 :               ico_l = coset(lxa, imax(lya - 1, 0), lza);
      95           0 :               jco_l = coset(lxb, (lyb + 1), lzb);
      96           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * tp->zeta[1] * pab;
      97             : 
      98           0 :               ico_l = coset(lxa, (lya + 1), lza);
      99           0 :               jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
     100           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) -= tp->zeta[0] * lyb * pab;
     101             : 
     102           0 :               ico_l = coset(lxa, (lya + 1), lza);
     103           0 :               jco_l = coset(lxb, (lyb + 1), lzb);
     104           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) +=
     105           0 :                   2.0 * tp->zeta[0] * tp->zeta[1] * pab;
     106             : 
     107             :               // z
     108             : 
     109           0 :               ico_l = coset(lxa, lya, imax(lza - 1, 0));
     110           0 :               jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
     111           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) += 0.5 * lza * lzb * pab;
     112             : 
     113           0 :               ico_l = coset(lxa, lya, imax(lza - 1, 0));
     114           0 :               jco_l = coset(lxb, lyb, (lzb + 1));
     115           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * tp->zeta[1] * pab;
     116             : 
     117           0 :               ico_l = coset(lxa, lya, (lza + 1));
     118           0 :               jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
     119           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) -= tp->zeta[0] * lzb * pab;
     120             : 
     121           0 :               ico_l = coset(lxa, lya, (lza + 1));
     122           0 :               jco_l = coset(lxb, lyb, (lzb + 1));
     123           0 :               idx2(tp->pab_prep[0], jco_l, ico_l) +=
     124           0 :                   2.0 * tp->zeta[0] * tp->zeta[1] * pab;
     125             :             }
     126             :           }
     127             :         }
     128             :       }
     129             :     }
     130             :   }
     131           0 : }
     132             : 
     133             : // *****************************************************************************
     134           0 : static void grid_prepare_pab_ADBmDAB(struct pab_computation_struct_ *const tp) {
     135             :   // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
     136             :   // is equivalent to mapping pab with
     137             :   //    pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b
     138             :   // ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) =
     139             :   //          pgf_a *(lbx pgf_{b-1x} - 2*tp->zeta[1]*pgf_{b+1x}) -
     140             :   //                   (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b
     141             : 
     142           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
     143           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
     144           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
     145           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
     146           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
     147           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
     148           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
     149           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
     150           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
     151           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
     152           0 :               const double pab = idx2(tp->pab[0], jco, ico);
     153           0 :               int ico_l, jco_l;
     154             : 
     155             :               // ! this element of pab results in 4 elements of pab_prep
     156           0 :               switch (tp->dir1) {
     157           0 :               case 'X': { // x
     158           0 :                 ico_l = coset(lxa, lya, lza);
     159           0 :                 jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
     160           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
     161             : 
     162             :                 // ico_l = coset(lxa, lya, lza);
     163           0 :                 jco_l = coset((lxb + 1), lyb, lzb);
     164           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
     165             : 
     166           0 :                 ico_l = coset(imax(lxa - 1, 0), lya, lza);
     167           0 :                 jco_l = coset(lxb, lyb, lzb);
     168           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
     169             : 
     170           0 :                 ico_l = coset(lxa + 1, lya, lza);
     171             :                 // jco_l = coset(lxb, lyb, lzb);
     172           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += 2.0 * tp->zeta[0] * pab;
     173           0 :               } break;
     174           0 :               case 'Y': { // y
     175           0 :                 ico_l = coset(lxa, lya, lza);
     176           0 :                 jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
     177           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
     178             : 
     179             :                 // ico_l = coset(lxa, lya, lza);
     180           0 :                 jco_l = coset(lxb, lyb + 1, lzb);
     181           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
     182             : 
     183           0 :                 ico_l = coset(lxa, imax(lya - 1, 0), lza);
     184           0 :                 jco_l = coset(lxb, lyb, lzb);
     185           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
     186             : 
     187           0 :                 ico_l = coset(lxa, lya + 1, lza);
     188             :                 // jco_l = coset(lxb, lyb, lzb);
     189           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += 2.0 * tp->zeta[0] * pab;
     190           0 :               } break;
     191           0 :               case 'Z': { // z
     192           0 :                 ico_l = coset(lxa, lya, lza);
     193           0 :                 jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
     194           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
     195             : 
     196             :                 // ico_l = coset(lxa, lya, lza);
     197           0 :                 jco_l = coset(lxb, lyb, lzb + 1);
     198           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
     199             : 
     200           0 :                 ico_l = coset(lxa, lya, imax(lza - 1, 0));
     201           0 :                 jco_l = coset(lxb, lyb, lzb);
     202           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
     203             : 
     204           0 :                 ico_l = coset(lxa, lya, lza + 1);
     205             :                 // jco_l = coset(lxb, lyb, lzb);
     206           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += 2.0 * tp->zeta[0] * pab;
     207           0 :               } break;
     208             :               default:
     209             :                 break;
     210             :               }
     211             :             }
     212             :           }
     213             :         }
     214             :       }
     215             :     }
     216             :   }
     217           0 : }
     218             : // *****************************************************************************
     219             : static void
     220           0 : grid_prepare_pab_ARDBmDARB(struct pab_computation_struct_ *const tp) {
     221             :   // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
     222             :   // is equivalent to mapping pab with
     223             :   // pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir}
     224             :   // pgf_b ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b )
     225             :   // =
     226             :   //                        pgf_a *(lbx pgf_{b-1x+1ir} -
     227             :   //                        2*tp->zeta[1]*pgf_{b+1x+1ir}) -
     228             :   //                       (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir}
     229             : 
     230           0 :   assert(1 <= tp->dir1 && tp->dir1 <= 3);
     231           0 :   assert(1 <= tp->dir2 && tp->dir2 <= 3);
     232             : 
     233           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
     234           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
     235           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
     236           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
     237           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
     238           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
     239           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
     240           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
     241           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
     242           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
     243           0 :               const double pab = idx2(tp->pab[0], jco, ico);
     244             : 
     245           0 :               int ico_l, jco_l;
     246             : 
     247             :               // this element of pab results in 4 elements of pab_prep
     248           0 :               switch (tp->dir1) {
     249           0 :               case 'X': {
     250           0 :                 switch (tp->dir2) {
     251           0 :                 case 'X': {
     252           0 :                   ico_l = coset(lxa, lya, lza);
     253           0 :                   jco_l = coset(lxb, lyb, lzb);
     254           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
     255             : 
     256           0 :                   ico_l = coset(lxa, lya, lza);
     257           0 :                   jco_l = coset((lxb + 2), lyb, lzb);
     258           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     259           0 :                       2.0 * tp->zeta[1] * pab;
     260             : 
     261           0 :                   ico_l = coset(imax(lxa - 1, 0), lya, lza);
     262           0 :                   jco_l = coset((lxb + 1), lyb, lzb);
     263           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
     264             : 
     265           0 :                   ico_l = coset((lxa + 1), lya, lza);
     266           0 :                   jco_l = coset((lxb + 1), lyb, lzb);
     267           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     268           0 :                       2.0 * tp->zeta[0] * pab;
     269           0 :                 } break;
     270           0 :                 case 'Y': {
     271           0 :                   ico_l = coset(lxa, lya, lza);
     272           0 :                   jco_l = coset(imax(lxb - 1, 0), (lyb + 1), lzb);
     273           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
     274             : 
     275           0 :                   ico_l = coset(lxa, lya, lza);
     276           0 :                   jco_l = coset((lxb + 1), (lyb + 1), lzb);
     277           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     278           0 :                       2.0 * tp->zeta[1] * pab;
     279             : 
     280           0 :                   ico_l = coset(imax(lxa - 1, 0), lya, lza);
     281           0 :                   jco_l = coset(lxb, (lyb + 1), lzb);
     282           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
     283             : 
     284           0 :                   ico_l = coset((lxa + 1), lya, lza);
     285           0 :                   jco_l = coset(lxb, (lyb + 1), lzb);
     286           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     287           0 :                       2.0 * tp->zeta[0] * pab;
     288           0 :                 } break;
     289           0 :                 case 'Z': {
     290           0 :                   ico_l = coset(lxa, lya, lza);
     291           0 :                   jco_l = coset(imax(lxb - 1, 0), lyb, (lzb + 1));
     292           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
     293             : 
     294           0 :                   ico_l = coset(lxa, lya, lza);
     295           0 :                   jco_l = coset((lxb + 1), lyb, (lzb + 1));
     296           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     297           0 :                       2.0 * tp->zeta[1] * pab;
     298             : 
     299           0 :                   ico_l = coset(imax(lxa - 1, 0), lya, lza);
     300           0 :                   jco_l = coset(lxb, lyb, (lzb + 1));
     301           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
     302             : 
     303           0 :                   ico_l = coset((lxa + 1), lya, lza);
     304           0 :                   jco_l = coset(lxb, lyb, (lzb + 1));
     305           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     306           0 :                       2.0 * tp->zeta[0] * pab;
     307           0 :                 } break;
     308             :                 default:
     309             :                   break;
     310             :                 }
     311             :               } break;
     312           0 :               case 'Y': {
     313           0 :                 switch (tp->dir2) {
     314           0 :                 case 'X': {
     315           0 :                   ico_l = coset(lxa, lya, lza);
     316           0 :                   jco_l = coset((lxb + 1), imax(lyb - 1, 0), lzb);
     317           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
     318             : 
     319           0 :                   ico_l = coset(lxa, lya, lza);
     320           0 :                   jco_l = coset((lxb + 1), (lyb + 1), lzb);
     321           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     322           0 :                       2.0 * tp->zeta[1] * pab;
     323             : 
     324           0 :                   ico_l = coset(lxa, imax(lya - 1, 0), lza);
     325           0 :                   jco_l = coset((lxb + 1), lyb, lzb);
     326           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
     327             : 
     328           0 :                   ico_l = coset(lxa, (lya + 1), lza);
     329           0 :                   jco_l = coset((lxb + 1), lyb, lzb);
     330           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     331           0 :                       2.0 * tp->zeta[0] * pab;
     332           0 :                 } break;
     333           0 :                 case 'Y': {
     334           0 :                   ico_l = coset(lxa, lya, lza);
     335           0 :                   jco_l = coset(lxb, lyb, lzb);
     336           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
     337             : 
     338           0 :                   ico_l = coset(lxa, lya, lza);
     339           0 :                   jco_l = coset(lxb, (lyb + 2), lzb);
     340           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     341           0 :                       2.0 * tp->zeta[1] * pab;
     342             : 
     343           0 :                   ico_l = coset(lxa, imax(lya - 1, 0), lza);
     344           0 :                   jco_l = coset(lxb, (lyb + 1), lzb);
     345           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
     346             : 
     347           0 :                   ico_l = coset(lxa, (lya + 1), lza);
     348           0 :                   jco_l = coset(lxb, (lyb + 1), lzb);
     349           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     350           0 :                       2.0 * tp->zeta[0] * pab;
     351           0 :                 } break;
     352           0 :                 case 'Z': {
     353           0 :                   ico_l = coset(lxa, lya, lza);
     354           0 :                   jco_l = coset(lxb, imax(lyb - 1, 0), (lzb + 1));
     355           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
     356             : 
     357           0 :                   ico_l = coset(lxa, lya, lza);
     358           0 :                   jco_l = coset(lxb, (lyb + 1), (lzb + 1));
     359           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     360           0 :                       2.0 * tp->zeta[1] * pab;
     361             : 
     362           0 :                   ico_l = coset(lxa, imax(lya - 1, 0), lza);
     363           0 :                   jco_l = coset(lxb, lyb, (lzb + 1));
     364           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
     365             : 
     366           0 :                   ico_l = coset(lxa, (lya + 1), lza);
     367           0 :                   jco_l = coset(lxb, lyb, (lzb + 1));
     368           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     369           0 :                       2.0 * tp->zeta[0] * pab;
     370           0 :                 } break;
     371             :                 default:
     372             :                   break;
     373             :                 }
     374             :               } break;
     375           0 :               case 'Z': {
     376           0 :                 switch (tp->dir2) {
     377           0 :                 case 'X': {
     378           0 :                   ico_l = coset(lxa, lya, lza);
     379           0 :                   jco_l = coset((lxb + 1), lyb, imax(lzb - 1, 0));
     380           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
     381             : 
     382           0 :                   ico_l = coset(lxa, lya, lza);
     383           0 :                   jco_l = coset((lxb + 1), lyb, (lzb + 1));
     384           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     385           0 :                       2.0 * tp->zeta[1] * pab;
     386             : 
     387           0 :                   ico_l = coset(lxa, lya, imax(lza - 1, 0));
     388           0 :                   jco_l = coset((lxb + 1), lyb, lzb);
     389           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
     390             : 
     391           0 :                   ico_l = coset(lxa, lya, (lza + 1));
     392           0 :                   jco_l = coset((lxb + 1), lyb, lzb);
     393           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     394           0 :                       2.0 * tp->zeta[0] * pab;
     395           0 :                 } break;
     396           0 :                 case 'Y': {
     397           0 :                   ico_l = coset(lxa, lya, lza);
     398           0 :                   jco_l = coset(lxb, (lyb + 1), imax(lzb - 1, 0));
     399           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
     400             : 
     401           0 :                   ico_l = coset(lxa, lya, lza);
     402           0 :                   jco_l = coset(lxb, (lyb + 1), (lzb + 1));
     403           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     404           0 :                       -2.0 * tp->zeta[1] * pab;
     405             : 
     406           0 :                   ico_l = coset(lxa, lya, imax(lza - 1, 0));
     407           0 :                   jco_l = coset(lxb, (lyb + 1), lzb);
     408           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
     409             : 
     410           0 :                   ico_l = coset(lxa, lya, (lza + 1));
     411           0 :                   jco_l = coset(lxb, (lyb + 1), lzb);
     412           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     413           0 :                       2.0 * tp->zeta[0] * pab;
     414           0 :                 } break;
     415           0 :                 case 'Z': {
     416           0 :                   ico_l = coset(lxa, lya, lza);
     417           0 :                   jco_l = coset(lxb, lyb, lzb);
     418           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
     419             : 
     420           0 :                   ico_l = coset(lxa, lya, lza);
     421           0 :                   jco_l = coset(lxb, lyb, (lzb + 2));
     422           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -=
     423           0 :                       2.0 * tp->zeta[1] * pab;
     424             : 
     425           0 :                   ico_l = coset(lxa, lya, imax(lza - 1, 0));
     426           0 :                   jco_l = coset(lxb, lyb, (lzb + 1));
     427           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
     428             : 
     429           0 :                   ico_l = coset(lxa, lya, (lza + 1));
     430           0 :                   jco_l = coset(lxb, lyb, (lzb + 1));
     431           0 :                   idx2(tp->pab_prep[0], jco_l, ico_l) +=
     432           0 :                       2.0 * tp->zeta[0] * pab;
     433           0 :                 } break;
     434             :                 default:
     435             :                   break;
     436             :                 }
     437             :               } break;
     438             :               default:
     439             :                 break;
     440             :               }
     441             :             }
     442             :           }
     443             :         }
     444             :       }
     445             :     }
     446             :   }
     447           0 : }
     448             : // *****************************************************************************
     449           0 : static void grid_prepare_pab_DABpADB(struct pab_computation_struct_ *const tp) {
     450             :   // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
     451             :   // is equivalent to mapping pab with
     452             :   //    pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b
     453             :   // ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) =
     454             :   //          pgf_a *(lbx pgf_{b-1x} + 2*tp->zeta[1]*pgf_{b+1x}) +
     455             :   //                   (lax pgf_{a-1x} + 2*zeta*pgf_{a+1x}) pgf_b
     456           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
     457           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
     458           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
     459           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
     460           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
     461           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
     462           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
     463           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
     464           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
     465           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
     466           0 :               const double pab = idx2(tp->pab[0], jco, ico);
     467             : 
     468           0 :               int ico_l, jco_l;
     469             : 
     470             :               // this element of pab results in 4 elements of pab_prep
     471             : 
     472           0 :               switch (tp->dir1) {
     473           0 :               case 'X': {
     474           0 :                 ico_l = coset(lxa, lya, lza);
     475           0 :                 jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
     476           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
     477             : 
     478           0 :                 ico_l = coset(lxa, lya, lza);
     479           0 :                 jco_l = coset((lxb + 1), lyb, lzb);
     480           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
     481             : 
     482           0 :                 ico_l = coset(imax(lxa - 1, 0), lya, lza);
     483           0 :                 jco_l = coset(lxb, lyb, lzb);
     484           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lxa * pab;
     485             : 
     486           0 :                 ico_l = coset((lxa + 1), lya, lza);
     487           0 :                 jco_l = coset(lxb, lyb, lzb);
     488           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[0] * pab;
     489           0 :               } break;
     490           0 :               case 'Y': { // y
     491           0 :                 ico_l = coset(lxa, lya, lza);
     492           0 :                 jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
     493           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
     494             : 
     495           0 :                 ico_l = coset(lxa, lya, lza);
     496           0 :                 jco_l = coset(lxb, (lyb + 1), lzb);
     497           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
     498             : 
     499           0 :                 ico_l = coset(lxa, imax(lya - 1, 0), lza);
     500           0 :                 jco_l = coset(lxb, lyb, lzb);
     501           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lya * pab;
     502             : 
     503           0 :                 ico_l = coset(lxa, (lya + 1), lza);
     504           0 :                 jco_l = coset(lxb, lyb, lzb);
     505           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[0] * pab;
     506           0 :               } break;
     507           0 :               case 'Z': { // z
     508           0 :                 ico_l = coset(lxa, lya, lza);
     509           0 :                 jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
     510           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
     511             : 
     512           0 :                 ico_l = coset(lxa, lya, lza);
     513           0 :                 jco_l = coset(lxb, lyb, (lzb + 1));
     514           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
     515             : 
     516           0 :                 ico_l = coset(lxa, lya, imax(lza - 1, 0));
     517           0 :                 jco_l = coset(lxb, lyb, lzb);
     518           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lza * pab;
     519             : 
     520           0 :                 ico_l = coset(lxa, lya, (lza + 1));
     521           0 :                 jco_l = coset(lxb, lyb, lzb);
     522           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[0] * pab;
     523           0 :                 break;
     524             :               }
     525             :               default:
     526             :                 break;
     527             :               }
     528             :             }
     529             :           }
     530             :         }
     531             :       }
     532             :     }
     533             :   }
     534           0 : }
     535             : // *****************************************************************************
     536           0 : static void grid_prepare_pab_Di(struct pab_computation_struct_ *const tp) {
     537             :   // create a new pab_local so that mapping pab_local with pgf_a pgf_b
     538             :   // is equivalent to mapping pab with
     539             :   //   d_{ider1} pgf_a d_{ider1} pgf_b
     540             :   // dx pgf_a dx pgf_b =
     541             :   //        (lax pgf_{a-1x})*(lbx pgf_{b-1x}) -
     542             :   //        2*tp->zeta[1]*lax*pgf_{a-1x}*pgf_{b+1x} -
     543             :   //         lbx pgf_{b-1x}*2*zeta*pgf_{a+1x}+ 4*zeta*zetab*pgf_{a+1x}pgf_{b+1x}
     544             : 
     545           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
     546           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
     547           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
     548           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
     549           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
     550           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
     551           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
     552           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
     553           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
     554           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
     555           0 :               const double pab = idx2(tp->pab[0], jco, ico);
     556             : 
     557           0 :               int ico_l, jco_l;
     558             :               // this element of pab results in 12 elements of pab_prep
     559             : 
     560           0 :               switch (tp->dir1) {
     561           0 :               case 'X': {
     562             :                 // x  (all safe if lxa = 0, as the spurious added terms have
     563             :                 // zero prefactor)
     564           0 :                 ico_l = coset(imax(lxa - 1, 0), lya, lza);
     565           0 :                 jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
     566           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lxa * lxb * pab;
     567             : 
     568           0 :                 ico_l = coset(imax(lxa - 1, 0), lya, lza);
     569           0 :                 jco_l = coset((lxb + 1), lyb, lzb);
     570           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -=
     571           0 :                     2.0 * lxa * tp->zeta[1] * pab;
     572             : 
     573           0 :                 ico_l = coset((lxa + 1), lya, lza);
     574           0 :                 jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
     575           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -=
     576           0 :                     2.0 * tp->zeta[0] * lxb * pab;
     577             : 
     578           0 :                 ico_l = coset((lxa + 1), lya, lza);
     579           0 :                 jco_l = coset((lxb + 1), lyb, lzb);
     580           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) +=
     581           0 :                     4.0 * tp->zeta[0] * tp->zeta[1] * pab;
     582           0 :               } break;
     583           0 :               case 'Y': {
     584             :                 // y
     585           0 :                 ico_l = coset(lxa, imax(lya - 1, 0), lza);
     586           0 :                 jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
     587           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lya * lyb * pab;
     588             : 
     589           0 :                 ico_l = coset(lxa, imax(lya - 1, 0), lza);
     590           0 :                 jco_l = coset(lxb, (lyb + 1), lzb);
     591           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -=
     592           0 :                     2.0 * lya * tp->zeta[1] * pab;
     593             : 
     594           0 :                 ico_l = coset(lxa, (lya + 1), lza);
     595           0 :                 jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
     596           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -=
     597           0 :                     2.0 * tp->zeta[0] * lyb * pab;
     598             : 
     599           0 :                 ico_l = coset(lxa, (lya + 1), lza);
     600           0 :                 jco_l = coset(lxb, (lyb + 1), lzb);
     601           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) +=
     602           0 :                     4.0 * tp->zeta[0] * tp->zeta[1] * pab;
     603           0 :               } break;
     604           0 :               case 'Z': {
     605             :                 // z
     606           0 :                 ico_l = coset(lxa, lya, imax(lza - 1, 0));
     607           0 :                 jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
     608           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) += lza * lzb * pab;
     609             : 
     610           0 :                 ico_l = coset(lxa, lya, imax(lza - 1, 0));
     611           0 :                 jco_l = coset(lxb, lyb, (lzb + 1));
     612           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -=
     613           0 :                     2.0 * lza * tp->zeta[1] * pab;
     614             : 
     615           0 :                 ico_l = coset(lxa, lya, (lza + 1));
     616           0 :                 jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
     617           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) -=
     618           0 :                     2.0 * tp->zeta[0] * lzb * pab;
     619             : 
     620           0 :                 ico_l = coset(lxa, lya, (lza + 1));
     621           0 :                 jco_l = coset(lxb, lyb, (lzb + 1));
     622           0 :                 idx2(tp->pab_prep[0], jco_l, ico_l) +=
     623           0 :                     4.0 * tp->zeta[0] * tp->zeta[1] * pab;
     624           0 :               } break;
     625             :               default:
     626             :                 break;
     627             :               }
     628             :             }
     629             :           }
     630             :         }
     631             :       }
     632             :     }
     633             :   }
     634           0 : }
     635             : 
     636             : // *****************************************************************************
     637           0 : static void oneterm_dijdij(const int idir, const double func_a, const int ico_l,
     638             :                            const int lx, const int ly, const int lz,
     639             :                            const double zet, tensor *const pab_prep) {
     640           0 :   int jco_l;
     641             : 
     642           0 :   switch (idir) {
     643           0 :   case 'X': {
     644           0 :     const int l1 = lx;
     645           0 :     const int l2 = ly;
     646           0 :     jco_l = coset(imax(lx - 1, 0), imax(ly - 1, 0), lz);
     647           0 :     idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
     648             : 
     649           0 :     jco_l = coset(lx + 1, imax(ly - 1, 0), lz);
     650           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
     651             : 
     652           0 :     jco_l = coset(imax(lx - 1, 0), ly + 1, lz);
     653           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
     654             : 
     655           0 :     jco_l = coset(lx + 1, ly + 1, lz);
     656           0 :     idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
     657           0 :   } break;
     658           0 :   case 'Y': {
     659           0 :     const int l1 = ly;
     660           0 :     const int l2 = lz;
     661           0 :     jco_l = coset(lx, imax(ly - 1, 0), imax(lz - 1, 0));
     662           0 :     idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
     663             : 
     664           0 :     jco_l = coset(lx, ly + 1, imax(lz - 1, 0));
     665           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
     666             : 
     667           0 :     jco_l = coset(lx, imax(ly - 1, 0), lz + 1);
     668           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
     669             : 
     670           0 :     jco_l = coset(lx, ly + 1, lz + 1);
     671           0 :     idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
     672           0 :   } break;
     673           0 :   case 'Z': {
     674           0 :     const int l1 = lz;
     675           0 :     const int l2 = lx;
     676           0 :     jco_l = coset(imax(lx - 1, 0), ly, imax(lz - 1, 0));
     677           0 :     idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
     678             : 
     679           0 :     jco_l = coset(imax(lx - 1, 0), ly, lz + 1);
     680           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
     681             : 
     682           0 :     jco_l = coset(lx + 1, ly, imax(lz - 1, 0));
     683           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
     684             : 
     685           0 :     jco_l = coset(lx + 1, ly, lz + 1);
     686           0 :     idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
     687           0 :   } break;
     688             :   default:
     689             :     break;
     690             :   }
     691           0 : }
     692             : // *****************************************************************************
     693           0 : static void grid_prepare_pab_DiDj(struct pab_computation_struct_ *const tp) {
     694             :   // create a new pab_local so that mapping pab_local with pgf_a pgf_b
     695             :   // is equivalent to mapping pab with
     696             :   //   d_{ider1} pgf_a d_{ider1} pgf_b
     697             : 
     698           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
     699           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
     700           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
     701           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
     702           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
     703           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
     704           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
     705           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
     706           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
     707           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
     708           0 :               const double pab = idx2(tp->pab[0], jco, ico);
     709             : 
     710           0 :               int ico_l;
     711           0 :               double func_a;
     712             : 
     713             :               // this element of pab results in 12 elements of pab_local
     714             : 
     715           0 :               if ((tp->dir1 == 'X' && tp->dir2 == 'Y') ||
     716             :                   (tp->dir1 == 'Y' && tp->dir2 == 'X')) {
     717             :                 // xy
     718           0 :                 ico_l = coset(imax(lxa - 1, 0), imax(lya - 1, 0), lza);
     719           0 :                 func_a = lxa * lya * pab;
     720           0 :                 oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     721             :                                tp->pab_prep);
     722             : 
     723           0 :                 ico_l = coset(lxa + 1, imax(lya - 1, 0), lza);
     724           0 :                 func_a = -2.0 * tp->zeta[0] * lya * pab;
     725           0 :                 oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     726             :                                tp->pab_prep);
     727             : 
     728           0 :                 ico_l = coset(imax(lxa - 1, 0), lya + 1, lza);
     729           0 :                 func_a = -2.0 * tp->zeta[0] * lxa * pab;
     730           0 :                 oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     731             :                                tp->pab_prep);
     732             : 
     733           0 :                 ico_l = coset(lxa + 1, lya + 1, lza);
     734           0 :                 func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
     735           0 :                 oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     736             :                                tp->pab_prep);
     737           0 :               } else if ((tp->dir1 == 'Y' && tp->dir2 == 'Z') ||
     738             :                          (tp->dir1 == 'Z' && tp->dir2 == 'Y')) {
     739             :                 // yz
     740           0 :                 ico_l = coset(lxa, imax(lya - 1, 0), imax(lza - 1, 0));
     741           0 :                 func_a = lya * lza * pab;
     742           0 :                 oneterm_dijdij('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     743             :                                tp->pab_prep);
     744             : 
     745           0 :                 ico_l = coset(lxa, lya + 1, imax(lza - 1, 0));
     746           0 :                 func_a = -2.0 * tp->zeta[0] * lza * pab;
     747           0 :                 oneterm_dijdij('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     748             :                                tp->pab_prep);
     749             : 
     750           0 :                 ico_l = coset(lxa, imax(lya - 1, 0), lza + 1);
     751           0 :                 func_a = -2.0 * tp->zeta[0] * lya * pab;
     752           0 :                 oneterm_dijdij('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     753             :                                tp->pab_prep);
     754             : 
     755           0 :                 ico_l = coset(lxa, lya + 1, lza + 1);
     756           0 :                 func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
     757           0 :                 oneterm_dijdij(2, func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     758             :                                tp->pab_prep);
     759           0 :               } else if ((tp->dir1 == 'Z' && tp->dir2 == 'X') ||
     760             :                          (tp->dir1 == 'X' && tp->dir2 == 'Z')) {
     761             :                 // zx
     762           0 :                 ico_l = coset(imax(lxa - 1, 0), lya, imax(lza - 1, 0));
     763           0 :                 func_a = lza * lxa * pab;
     764           0 :                 oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     765             :                                tp->pab_prep);
     766             : 
     767           0 :                 ico_l = coset(imax(lxa - 1, 0), lya, lza + 1);
     768           0 :                 func_a = -2.0 * tp->zeta[0] * lxa * pab;
     769           0 :                 oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     770             :                                tp->pab_prep);
     771             : 
     772           0 :                 ico_l = coset(lxa + 1, lya, imax(lza - 1, 0));
     773           0 :                 func_a = -2.0 * tp->zeta[0] * lza * pab;
     774           0 :                 oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     775             :                                tp->pab_prep);
     776             : 
     777           0 :                 ico_l = coset(lxa + 1, lya, lza + 1);
     778           0 :                 func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
     779           0 :                 oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     780             :                                tp->pab_prep);
     781             :               }
     782             :             }
     783             :           }
     784             :         }
     785             :       }
     786             :     }
     787             :   }
     788           0 : }
     789             : 
     790             : // *****************************************************************************
     791           0 : static void oneterm_diidii(const int idir, const double func_a, const int ico_l,
     792             :                            const int lx, const int ly, const int lz,
     793             :                            const double zet, tensor *const pab_prep) {
     794           0 :   int jco_l;
     795             : 
     796           0 :   switch (idir) {
     797           0 :   case 'X': {
     798           0 :     const int l1 = lx;
     799           0 :     jco_l = coset(imax(lx - 2, 0), ly, lz);
     800           0 :     idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
     801             : 
     802           0 :     jco_l = coset(lx, ly, lz);
     803           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
     804             : 
     805           0 :     jco_l = coset(lx + 2, ly, lz);
     806           0 :     idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
     807           0 :   } break;
     808           0 :   case 'Y': {
     809           0 :     const int l1 = ly;
     810           0 :     jco_l = coset(lx, imax(ly - 2, 0), lz);
     811           0 :     idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
     812             : 
     813           0 :     jco_l = coset(lx, ly, lz);
     814           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
     815             : 
     816           0 :     jco_l = coset(lx, ly + 2, lz);
     817           0 :     idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
     818           0 :   } break;
     819           0 :   case 'Z': {
     820           0 :     const int l1 = lz;
     821           0 :     jco_l = coset(lx, ly, imax(lz - 2, 0));
     822           0 :     idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
     823             : 
     824           0 :     jco_l = coset(lx, ly, lz);
     825           0 :     idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
     826             : 
     827           0 :     jco_l = coset(lx, ly, lz + 2);
     828           0 :     idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
     829           0 :   } break;
     830             :   default:
     831           0 :     printf("Wrong value for ider: should be 1, 2, or 3\n");
     832           0 :     abort();
     833           0 :     break;
     834             :   }
     835           0 : }
     836             : 
     837             : // *****************************************************************************
     838           0 : static void grid_prepare_pab_Di2(struct pab_computation_struct_ *const tp) {
     839             :   // create a new pab_local so that mapping pab_local with pgf_a pgf_b
     840             :   // is equivalent to mapping pab with
     841             :   //   dd_{ider1} pgf_a dd_{ider1} pgf_b
     842             : 
     843           0 :   for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
     844           0 :     for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
     845           0 :       for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
     846           0 :         for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
     847           0 :           for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
     848           0 :                lza <= tp->lmax[0] - lxa - lya; lza++) {
     849           0 :             for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
     850           0 :                  lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
     851           0 :               const int ico = tp->offset[0] + coset(lxa, lya, lza);
     852           0 :               const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
     853           0 :               const double pab = idx2(tp->pab[0], jco, ico);
     854             : 
     855           0 :               int ico_l;
     856           0 :               double func_a;
     857             : 
     858             :               // this element of pab results in  9 elements of pab_local
     859           0 :               switch (tp->dir1) {
     860           0 :               case 'X': {
     861             :                 // x
     862           0 :                 ico_l = coset(imax(lxa - 2, 0), lya, lza);
     863           0 :                 func_a = lxa * (lxa - 1) * pab;
     864           0 :                 oneterm_diidii('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     865             :                                tp->pab_prep);
     866             : 
     867           0 :                 ico_l = coset(lxa, lya, lza);
     868           0 :                 func_a = -2.0 * tp->zeta[0] * (2 * lxa + 1) * pab;
     869           0 :                 oneterm_diidii('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     870             :                                tp->pab_prep);
     871             : 
     872           0 :                 ico_l = coset(lxa + 2, lya, lza);
     873           0 :                 func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
     874           0 :                 oneterm_diidii('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     875             :                                tp->pab_prep);
     876           0 :               } break;
     877           0 :               case 'Y': {
     878             :                 // y
     879           0 :                 ico_l = coset(lxa, imax(lya - 2, 0), lza);
     880           0 :                 func_a = lya * (lya - 1) * pab;
     881           0 :                 oneterm_diidii('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     882             :                                tp->pab_prep);
     883             : 
     884           0 :                 ico_l = coset(lxa, lya, lza);
     885           0 :                 func_a = -2.0 * tp->zeta[0] * (2 * lya + 1) * pab;
     886           0 :                 oneterm_diidii('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     887             :                                tp->pab_prep);
     888             : 
     889           0 :                 ico_l = coset(lxa, lya + 2, lza);
     890           0 :                 func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
     891           0 :                 oneterm_diidii('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     892             :                                tp->pab_prep);
     893           0 :               } break;
     894           0 :               case 'Z': {
     895             :                 // z
     896           0 :                 ico_l = coset(lxa, lya, imax(lza - 2, 0));
     897           0 :                 func_a = lza * (lza - 1) * pab;
     898           0 :                 oneterm_diidii('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     899             :                                tp->pab_prep);
     900             : 
     901           0 :                 ico_l = coset(lxa, lya, lza);
     902           0 :                 func_a = -2.0 * tp->zeta[0] * (2 * lza + 1) * pab;
     903           0 :                 oneterm_diidii('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     904             :                                tp->pab_prep);
     905             : 
     906           0 :                 ico_l = coset(lxa, lya, lza + 2);
     907           0 :                 func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
     908           0 :                 oneterm_diidii('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
     909             :                                tp->pab_prep);
     910           0 :               } break;
     911             :               default:
     912           0 :                 printf("Wrong value for ider: should be 'X', 'Y' or 'Z'.\n");
     913           0 :                 abort();
     914           0 :                 break;
     915             :               }
     916             :             }
     917             :           }
     918             :         }
     919             :       }
     920             :     }
     921             :   }
     922           0 : }
     923             : 
     924             : // *****************************************************************************
     925         632 : void grid_prepare_get_ldiffs_dgemm(const enum grid_func func,
     926             :                                    int *const lmin_diff, int *const lmax_diff) {
     927         632 :   switch (func) {
     928         632 :   case GRID_FUNC_AB:
     929         632 :     lmax_diff[0] = 0;
     930         632 :     lmin_diff[0] = 0;
     931         632 :     lmax_diff[1] = 0;
     932         632 :     lmin_diff[1] = 0;
     933         632 :     break;
     934           0 :   case GRID_FUNC_DADB:
     935             :   case GRID_FUNC_ADBmDAB_X:
     936             :   case GRID_FUNC_ADBmDAB_Y:
     937             :   case GRID_FUNC_ADBmDAB_Z:
     938             :   case GRID_FUNC_DABpADB_X:
     939             :   case GRID_FUNC_DABpADB_Y:
     940             :   case GRID_FUNC_DABpADB_Z:
     941           0 :     lmax_diff[0] = 1;
     942           0 :     lmin_diff[0] = -1;
     943           0 :     lmax_diff[1] = 1;
     944           0 :     lmin_diff[1] = -1;
     945           0 :     break;
     946           0 :   case GRID_FUNC_ARDBmDARB_XX:
     947             :   case GRID_FUNC_ARDBmDARB_XY:
     948             :   case GRID_FUNC_ARDBmDARB_XZ:
     949             :   case GRID_FUNC_ARDBmDARB_YX:
     950             :   case GRID_FUNC_ARDBmDARB_YY:
     951             :   case GRID_FUNC_ARDBmDARB_YZ:
     952             :   case GRID_FUNC_ARDBmDARB_ZX:
     953             :   case GRID_FUNC_ARDBmDARB_ZY:
     954             :   case GRID_FUNC_ARDBmDARB_ZZ:
     955           0 :     lmax_diff[0] = 1; // TODO: mistake???, then we could merge la and lb.
     956           0 :     lmin_diff[0] = -1;
     957           0 :     lmax_diff[1] = 2;
     958           0 :     lmin_diff[1] = -1;
     959           0 :     break;
     960           0 :   case GRID_FUNC_DX:
     961             :   case GRID_FUNC_DY:
     962             :   case GRID_FUNC_DZ:
     963           0 :     lmax_diff[0] = 1;
     964           0 :     lmin_diff[0] = -1;
     965           0 :     lmax_diff[1] = 1;
     966           0 :     lmin_diff[1] = -1;
     967           0 :     break;
     968           0 :   case GRID_FUNC_DXDY:
     969             :   case GRID_FUNC_DYDZ:
     970             :   case GRID_FUNC_DZDX:
     971             :   case GRID_FUNC_DXDX:
     972             :   case GRID_FUNC_DYDY:
     973             :   case GRID_FUNC_DZDZ:
     974           0 :     lmax_diff[0] = 2;
     975           0 :     lmin_diff[0] = -2;
     976           0 :     lmax_diff[1] = 2;
     977           0 :     lmin_diff[1] = -2;
     978           0 :     break;
     979             :   default:
     980           0 :     printf("Unkown ga-gb function");
     981           0 :     abort();
     982             :   }
     983         632 : }
     984             : 
     985             : // *****************************************************************************
     986       44389 : void grid_prepare_pab_dgemm(const enum grid_func func, const int *const offset,
     987             :                             const int *const lmin, const int *const lmax,
     988             :                             const double *const zeta, tensor *const pab,
     989             :                             tensor *const pab_prep) {
     990       44389 :   struct pab_computation_struct_ tmp;
     991             : 
     992       44389 :   tmp.offset[0] = offset[0];
     993       44389 :   tmp.offset[1] = offset[1];
     994             : 
     995       44389 :   tmp.lmin[0] = lmin[0];
     996       44389 :   tmp.lmin[1] = lmin[1];
     997             : 
     998       44389 :   tmp.lmax[0] = lmax[0];
     999       44389 :   tmp.lmax[1] = lmax[1];
    1000             : 
    1001       44389 :   tmp.pab = pab;
    1002       44389 :   tmp.pab_prep = pab_prep;
    1003             : 
    1004       44389 :   tmp.zeta[0] = zeta[0];
    1005       44389 :   tmp.zeta[1] = zeta[1];
    1006       44389 :   memset(pab_prep->data, 0, pab_prep->alloc_size_ * sizeof(double));
    1007             : 
    1008       44389 :   switch (func) {
    1009       44389 :   case GRID_FUNC_AB:
    1010       44389 :     grid_prepare_pab_AB(&tmp);
    1011       44389 :     break;
    1012           0 :   case GRID_FUNC_DADB:
    1013           0 :     grid_prepare_pab_DADB(&tmp);
    1014           0 :     break;
    1015           0 :   case GRID_FUNC_ADBmDAB_X:
    1016           0 :     tmp.dir1 = 'X';
    1017           0 :     grid_prepare_pab_ADBmDAB(&tmp);
    1018           0 :     break;
    1019           0 :   case GRID_FUNC_ADBmDAB_Y:
    1020           0 :     tmp.dir1 = 'Y';
    1021           0 :     grid_prepare_pab_ADBmDAB(&tmp);
    1022           0 :     break;
    1023           0 :   case GRID_FUNC_ADBmDAB_Z:
    1024           0 :     tmp.dir1 = 'Z';
    1025           0 :     grid_prepare_pab_ADBmDAB(&tmp);
    1026           0 :     break;
    1027           0 :   case GRID_FUNC_ARDBmDARB_XX:
    1028           0 :     tmp.dir1 = 'X';
    1029           0 :     tmp.dir2 = 'X';
    1030           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1031           0 :     break;
    1032           0 :   case GRID_FUNC_ARDBmDARB_XY:
    1033           0 :     tmp.dir1 = 'X';
    1034           0 :     tmp.dir2 = 'Y';
    1035           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1036           0 :     break;
    1037           0 :   case GRID_FUNC_ARDBmDARB_XZ:
    1038           0 :     tmp.dir1 = 'X';
    1039           0 :     tmp.dir2 = 'Z';
    1040           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1041           0 :     break;
    1042           0 :   case GRID_FUNC_ARDBmDARB_YX:
    1043           0 :     tmp.dir1 = 'Y';
    1044           0 :     tmp.dir2 = 'X';
    1045           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1046           0 :     break;
    1047           0 :   case GRID_FUNC_ARDBmDARB_YY:
    1048           0 :     tmp.dir1 = 'Y';
    1049           0 :     tmp.dir2 = 'Y';
    1050           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1051           0 :     break;
    1052           0 :   case GRID_FUNC_ARDBmDARB_YZ:
    1053           0 :     tmp.dir1 = 'Y';
    1054           0 :     tmp.dir2 = 'Z';
    1055           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1056           0 :     break;
    1057           0 :   case GRID_FUNC_ARDBmDARB_ZX:
    1058           0 :     tmp.dir1 = 'Z';
    1059           0 :     tmp.dir2 = 'X';
    1060           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1061           0 :     break;
    1062           0 :   case GRID_FUNC_ARDBmDARB_ZY:
    1063           0 :     tmp.dir1 = 'Z';
    1064           0 :     tmp.dir2 = 'Y';
    1065           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1066           0 :     break;
    1067           0 :   case GRID_FUNC_ARDBmDARB_ZZ:
    1068           0 :     tmp.dir1 = 'Z';
    1069           0 :     tmp.dir2 = 'Z';
    1070           0 :     grid_prepare_pab_ARDBmDARB(&tmp);
    1071           0 :     break;
    1072           0 :   case GRID_FUNC_DABpADB_X:
    1073           0 :     tmp.dir1 = 'X';
    1074           0 :     grid_prepare_pab_DABpADB(&tmp);
    1075           0 :     break;
    1076           0 :   case GRID_FUNC_DABpADB_Y:
    1077           0 :     tmp.dir1 = 'Y';
    1078           0 :     grid_prepare_pab_DABpADB(&tmp);
    1079           0 :     break;
    1080           0 :   case GRID_FUNC_DABpADB_Z:
    1081           0 :     tmp.dir1 = 'Z';
    1082           0 :     grid_prepare_pab_DABpADB(&tmp);
    1083           0 :     break;
    1084           0 :   case GRID_FUNC_DX:
    1085           0 :     tmp.dir1 = 'X';
    1086           0 :     grid_prepare_pab_Di(&tmp);
    1087           0 :     break;
    1088           0 :   case GRID_FUNC_DY:
    1089           0 :     tmp.dir1 = 'Y';
    1090           0 :     grid_prepare_pab_Di(&tmp);
    1091           0 :     break;
    1092           0 :   case GRID_FUNC_DZ:
    1093           0 :     tmp.dir1 = 'Z';
    1094           0 :     grid_prepare_pab_Di(&tmp);
    1095           0 :     break;
    1096           0 :   case GRID_FUNC_DXDY:
    1097           0 :     tmp.dir1 = 'X';
    1098           0 :     tmp.dir2 = 'Y';
    1099           0 :     grid_prepare_pab_DiDj(&tmp);
    1100           0 :     break;
    1101           0 :   case GRID_FUNC_DYDZ:
    1102           0 :     tmp.dir1 = 'Y';
    1103           0 :     tmp.dir2 = 'Z';
    1104           0 :     grid_prepare_pab_DiDj(&tmp);
    1105           0 :     break;
    1106           0 :   case GRID_FUNC_DZDX:
    1107           0 :     tmp.dir1 = 'Z';
    1108           0 :     tmp.dir2 = 'X';
    1109           0 :     grid_prepare_pab_DiDj(&tmp);
    1110           0 :     break;
    1111           0 :   case GRID_FUNC_DXDX:
    1112           0 :     tmp.dir1 = 'X';
    1113           0 :     grid_prepare_pab_Di2(&tmp);
    1114           0 :     break;
    1115           0 :   case GRID_FUNC_DYDY:
    1116           0 :     tmp.dir1 = 'Y';
    1117           0 :     grid_prepare_pab_Di2(&tmp);
    1118           0 :     break;
    1119           0 :   case GRID_FUNC_DZDZ:
    1120           0 :     tmp.dir1 = 'Z';
    1121           0 :     grid_prepare_pab_Di2(&tmp);
    1122           0 :     break;
    1123             :   default:
    1124           0 :     printf("Unkown ga-gb function");
    1125           0 :     abort();
    1126             :   }
    1127       44389 : }

Generated by: LCOV version 1.15