ode/src/threaded_solver_ldlt.h

   1 /*************************************************************************
   2  *                                                                       *
   3  * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith.       *
   4  * All rights reserved.  Email: russ@q12.org   Web: www.q12.org          *
   5  *                                                                       *
   6  * This library is free software; you can redistribute it and/or         *
   7  * modify it under the terms of EITHER:                                  *
   8  *   (1) The GNU Lesser General Public License as published by the Free  *
   9  *       Software Foundation; either version 2.1 of the License, or (at  *
  10  *       your option) any later version. The text of the GNU Lesser      *
  11  *       General Public License is included with this library in the     *
  12  *       file LICENSE.TXT.                                               *
  13  *   (2) The BSD-style license that is included with this library in     *
  14  *       the file LICENSE-BSD.TXT.                                       *
  15  *                                                                       *
  16  * This library is distributed in the hope that it will be useful,       *
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files    *
  19  * LICENSE.TXT and LICENSE-BSD.TXT for more details.                     *
  20  *                                                                       *
  21  *************************************************************************/
  22
  23 /*
  24  * Equation System Threaded Solver
  25  * Copyright (c) 2017-2024 Oleh Derevenko, odar@eleks.com (change all "a" to "e")
  26  */
  27
  28
  29
  30 #ifndef _ODE_THREADED_SOLVER_LDLT_H_
  31 #define _ODE_THREADED_SOLVER_LDLT_H_
  32
  33
  34 #include "coop_matrix_types.h"
  35 #include <ode/threading.h>
  36
  37
  38 class dxThreadingBase;
  39 class dxResourceRequirementDescriptor;
  40 class dxRequiredResourceContainer;
  41
  42
  43 class ThreadedEquationSolverLDLT
  44 {
  45 public:
  46     static void estimateCooperativeFactoringLDLTResourceRequirements(dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
  47         unsigned allowedThreadCount, unsigned rowCount);
  48     static void cooperativelyFactorLDLT(dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
  49         dReal *A, dReal *d, unsigned rowCount, unsigned rowSkip);
  50
  51     static void estimateCooperativeSolvingL1StraightResourceRequirements(dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
  52         unsigned allowedThreadCount, unsigned rowCount);
  53     static void cooperativelySolveL1Straight(dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
  54         const dReal *L, dReal *b, unsigned rowCount, unsigned rowSkip);
  55
  56     static void estimateCooperativeSolvingL1TransposedResourceRequirements(dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
  57         unsigned allowedThreadCount, unsigned rowCount);
  58     static void cooperativelySolveL1Transposed(dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
  59         const dReal *L, dReal *b, unsigned rowCount, unsigned rowSkip);
  60
  61     static void estimateCooperativeScalingVectorResourceRequirements(dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
  62         unsigned allowedThreadCount, unsigned elementCount);
  63     static void cooperativelyScaleVector(dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
  64         dReal *vectorData, const dReal *scaleData, unsigned elementCount);
  65
  66     static void estimateCooperativeSolvingLDLTResourceRequirements(dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
  67         unsigned allowedThreadCount, unsigned rowCount);
  68     static void cooperativelySolveLDLT(dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
  69         const dReal *L, const dReal *d, dReal *b, unsigned rowCount, unsigned rowSkip);
  70
  71 public:
  72     enum
  73     {
  74         ALLOCATION_DEFAULT_ALIGNMENT = COOP_THREAD_DATA_ALIGNMENT_SIZE,
  75     };
  76
  77 private:
  78     struct FactorizationSolveL1StripeCellContext;
  79     struct FactorizationFactorizeL1StripeThreadContext;
  80
  81     enum
  82     {
  83         FLDLT_D_STRIDE          = 1,
  84         FLDLT_COOPERATIVE_BLOCK_COUNT_MINIMUM = 5,
  85
  86         FSL1S_BLOCK_SIZE        = 2,
  87
  88         FSL1S_REGULAR_B_ROWS    = FSL1S_BLOCK_SIZE,
  89         FSL1S_FINAL_B_ROWS      = 1,
  90
  91         FFL1S_REGULAR_A_ROWS    = FSL1S_BLOCK_SIZE,
  92         FFL1S_FINAL_A_ROWS      = 1,
  93         FFL1S_REGULAR_BLOCK_SIZE = 16,  // A suitable by magnitude number being a power of 2 and (naturally) not being divisible by 6
  94         FFL1S_FINAL_BLOCK_SIZE  = 32, // A suitable by magnitude number being a power of 2 and (naturally) not being divisible by 6
  95     };
  96
  97     static unsigned restrictFactoringLDLTAllowedThreadCount(
  98         dxThreadingBase *threading, unsigned allowedThreadCount, unsigned rowCount);
  99     static void doEstimateCooperativeFactoringLDLTResourceRequirementsValidated(
 100         dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
 101         unsigned allowedThreadCount, unsigned rowCount);
 102     static void doCooperativelyFactorLDLTValidated(
 103         dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
 104         dReal *A, dReal *d, unsigned rowCount, unsigned rowSkip);
 105
 106
 107     static unsigned deriveSolvingL1StripeBlockCount(unsigned rowCount, unsigned blockStep)
 108     {
 109         return (rowCount + (blockStep - 1)) / blockStep;
 110     }
 111
 112     struct FactorizationSolvingL1StripeMemoryEstimates
 113     {
 114         void assignData(sizeint descriptorSizeRequired, sizeint contextSizeRequired)
 115         {
 116             m_descriptorSizeRequired = descriptorSizeRequired;
 117             m_contextSizeRequired = contextSizeRequired;
 118         }
 119
 120         sizeint  m_descriptorSizeRequired;
 121         sizeint  m_contextSizeRequired;
 122     };
 123
 124     static unsigned deriveSolvingL1StripeThreadCount(unsigned blockCount, unsigned allowedThreadCount)
 125     {
 126         dIASSERT(allowedThreadCount >= 1);
 127
 128         unsigned maximumCount = blockCount / 2;
 129         return maximumCount >= allowedThreadCount ? allowedThreadCount : dMACRO_MAX(maximumCount, 1U);
 130     }
 131
 132     static sizeint estimateCooperativelySolvingL1Stripe_XMemoryRequirement(unsigned blockCount,
 133         FactorizationSolvingL1StripeMemoryEstimates &ref_memoryEstimates)
 134     {
 135         sizeint descriptorSizeRequired = dOVERALIGNED_SIZE(sizeof(cellindexint) * blockCount, COOP_THREAD_DATA_ALIGNMENT_SIZE);
 136         sizeint contextSizeRequired = dOVERALIGNED_SIZE(sizeof(FactorizationSolveL1StripeCellContext) * (CCI__MAX + 1) * blockCount, COOP_THREAD_DATA_ALIGNMENT_SIZE);
 137         ref_memoryEstimates.assignData(descriptorSizeRequired, contextSizeRequired);
 138
 139         sizeint totalSizeRequired = descriptorSizeRequired + contextSizeRequired;
 140         return totalSizeRequired;
 141     }
 142
 143     static void *markCooperativelySolvingL1Stripe_XMemoryStructuresOut(void *buffer,
 144         const FactorizationSolvingL1StripeMemoryEstimates &memoryEstimates,
 145         cellindexint *&out_blockProgressDescriptors, FactorizationSolveL1StripeCellContext *&out_cellContexts)
 146     {
 147         void *currentLocation = buffer;
 148
 149         out_blockProgressDescriptors = (cellindexint *)currentLocation; currentLocation = (uint8 *)currentLocation + memoryEstimates.m_descriptorSizeRequired;
 150         out_cellContexts = (FactorizationSolveL1StripeCellContext *)currentLocation; currentLocation = (uint8 *)currentLocation + memoryEstimates.m_contextSizeRequired;
 151
 152         return currentLocation;
 153     }
 154
 155     static void initializeCooperativelySolvingL1Stripe_XMemoryStructures(unsigned blockCount,
 156         atomicord32 &out_blockCompletionProgress, cellindexint *blockProgressDescriptors, FactorizationSolveL1StripeCellContext *dUNUSED(cellContexts))
 157     {
 158         out_blockCompletionProgress = 0;
 159         memset(blockProgressDescriptors, 0, blockCount * sizeof(*blockProgressDescriptors));
 160     }
 161
 162     template<unsigned int block_step, unsigned int b_rows>
 163     static void participateSolvingL1Stripe_X(const dReal *L, dReal *B, unsigned blockCount, unsigned rowSkip,
 164         volatile atomicord32 &refBlockCompletionProgress/*=0*/, volatile cellindexint *blockProgressDescriptors/*=[blockCount]*/,
 165         FactorizationSolveL1StripeCellContext *cellContexts/*=[CCI__MAX x blockCount] + [blockCount]*/, unsigned ownThreadIndex);
 166
 167     static unsigned deriveScalingAndFactorizingL1StripeBlockCountFromSolvingBlockIndex(unsigned solvingBlockIndex, unsigned solvingBlockStep, unsigned blockARows)
 168     {
 169         unsigned factorizingBlockSize = deriveScalingAndFactorizingL1StripeBlockSize(blockARows);
 170         return deriveScalingAndFactorizingL1StripeBlockCountFromFactorizationRow(solvingBlockIndex * solvingBlockStep, factorizingBlockSize);
 171     }
 172
 173     static unsigned deriveScalingAndFactorizingL1StripeBlockCountFromFactorizationRow(unsigned factorizationRowIndex, unsigned factorizationBlockSize)
 174     {
 175         return (factorizationRowIndex + (factorizationBlockSize - 1)) / factorizationBlockSize;
 176     }
 177
 178     static unsigned deriveScalingAndFactorizingL1StripeBlockSize(unsigned blockARows)
 179     {
 180         unsigned result = blockARows != 1 ? FFL1S_REGULAR_BLOCK_SIZE : FFL1S_FINAL_BLOCK_SIZE;
 181         dIASSERT(blockARows >= 1 && blockARows <= 2);
 182
 183         return result;
 184     }
 185
 186
 187     static unsigned deriveScalingAndFactorizingL1StripeThreadCount(unsigned blockCount, unsigned allowedThreadCount)
 188     {
 189         dIASSERT(blockCount != 0);
 190         dIASSERT(allowedThreadCount >= 1);
 191
 192         return dMACRO_MIN(blockCount, allowedThreadCount);
 193     }
 194
 195     struct FactorizationFactorizeL1StripeContext;
 196
 197     struct FactorizationScalingAndFactorizingL1StripeMemoryEstimates
 198     {
 199         void assignData(sizeint contextSizeRequired)
 200         {
 201             m_contextSizeRequired = contextSizeRequired;
 202         }
 203
 204         sizeint  m_contextSizeRequired;
 205     };
 206
 207     static sizeint estimateCooperativelyScalingAndFactorizingL1Stripe_XMemoryRequirement(unsigned factorizingMaximumThreads,
 208         FactorizationScalingAndFactorizingL1StripeMemoryEstimates &ref_memoryEstimates)
 209     {
 210         dIASSERT(factorizingMaximumThreads != 0);
 211
 212         sizeint contextSizeRequired = dOVERALIGNED_SIZE(sizeof(FactorizationFactorizeL1StripeContext) + sizeof(FactorizationFactorizeL1StripeThreadContext) * (factorizingMaximumThreads - 1), COOP_THREAD_DATA_ALIGNMENT_SIZE);
 213         ref_memoryEstimates.assignData(contextSizeRequired);
 214
 215         sizeint totalSizeRequired = contextSizeRequired;
 216         return totalSizeRequired;
 217     }
 218
 219     static void *markCooperativelyScalingAndFactorizingL1Stripe_XMemoryStructuresOut(void *buffer,
 220         const FactorizationScalingAndFactorizingL1StripeMemoryEstimates &memoryEstimates, FactorizationFactorizeL1StripeContext *&out_factorizationContext)
 221     {
 222         void *currentLocation = buffer;
 223
 224         out_factorizationContext = (FactorizationFactorizeL1StripeContext *)currentLocation; currentLocation = (uint8 *)currentLocation + memoryEstimates.m_contextSizeRequired;
 225
 226         return currentLocation;
 227     }
 228
 229     static void initializeCooperativelyScalingAndFactorizingL1Stripe_XMemoryStructures(
 230         FactorizationFactorizeL1StripeContext *factorizationContext, unsigned threadCount)
 231     {
 232         factorizationContext->initialize(threadCount);
 233     }
 234
 235
 236     template<unsigned int a_rows, unsigned int d_stride>
 237     static void participateScalingAndFactorizingL1Stripe_X(dReal *ARow, dReal *d, unsigned factorizationRow, unsigned rowSkip,
 238         FactorizationFactorizeL1StripeContext *factorizationContext, unsigned ownThreadIndex);
 239
 240 private:
 241     struct FactorLDLTWorkerContext
 242     {
 243         FactorLDLTWorkerContext(dxThreadingBase *threading, unsigned allowedThreadCount,
 244             dReal *A, dReal *d, unsigned totalBlockCount, unsigned rowCount, unsigned rowSkip,
 245             atomicord32 &ref_solvingBlockCompletionProgress, cellindexint *solvingBlockProgressDescriptors,
 246             FactorizationSolveL1StripeCellContext *solvingCellContexts,
 247             FactorizationFactorizeL1StripeContext *factorizingFactorizationContext,
 248             dCallReleaseeID calculationFinishReleasee):
 249             m_threading(threading),
 250             m_allowedThreadCount(allowedThreadCount),
 251             m_A(A),
 252             m_ARow(A),
 253             m_d(d),
 254             m_solvingBlockIndex(0),
 255             m_totalBlockCount(totalBlockCount),
 256             m_rowCount(rowCount),
 257             m_rowSkip(rowSkip),
 258             m_refSolvingBlockCompletionProgress(ref_solvingBlockCompletionProgress),
 259             m_solvingBlockProgressDescriptors(solvingBlockProgressDescriptors),
 260             m_solvingCellContexts(solvingCellContexts),
 261             m_factorizingFactorizationContext(factorizingFactorizationContext),
 262             m_calculationFinishReleasee(calculationFinishReleasee)
 263         {
 264         }
 265
 266         void incrementForNextBlock()
 267         {
 268             const unsigned blockStep = FSL1S_BLOCK_SIZE;
 269
 270             m_ARow += blockStep * m_rowSkip;
 271             m_solvingBlockIndex += 1;
 272         }
 273
 274         dxThreadingBase             *m_threading;
 275         unsigned                    m_allowedThreadCount;
 276         dReal                       *m_A;
 277         dReal                       *m_ARow;
 278         dReal                       *m_d;
 279         unsigned                    m_solvingBlockIndex;
 280         unsigned                    m_totalBlockCount;
 281         unsigned                    m_rowCount;
 282         unsigned                    m_rowSkip;
 283         atomicord32                 &m_refSolvingBlockCompletionProgress;
 284         cellindexint                *m_solvingBlockProgressDescriptors;
 285         FactorizationSolveL1StripeCellContext *m_solvingCellContexts;
 286         FactorizationFactorizeL1StripeContext *m_factorizingFactorizationContext;
 287         dCallReleaseeID             m_calculationFinishReleasee;
 288     };
 289
 290     static int factotLDLT_solvingComplete_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 291     static void factotLDLT_solvingComplete(FactorLDLTWorkerContext &ref_context, unsigned ownThreadIndex);
 292
 293     static int factotLDLT_solvingCompleteSync_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 294     static void factotLDLT_solvingCompleteSync(FactorLDLTWorkerContext &ref_workerContext);
 295
 296     static int factotLDLT_scalingAndFactorizingComplete_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 297     static void factotLDLT_scalingAndFactorizingComplete(FactorLDLTWorkerContext &ref_workerContext, unsigned ownThreadIndex);
 298
 299     static int factotLDLT_scalingAndFactorizingCompleteSync_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 300     static void factotLDLT_scalingAndFactorizingCompleteSync(FactorLDLTWorkerContext &ref_workerContext);
 301
 302     static int factotLDLT_solvingFinal_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 303     static void factotLDLT_solvingFinal(FactorLDLTWorkerContext &ref_context, unsigned ownThreadIndex);
 304
 305     static int factotLDLT_solvingFinalSync_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 306     static void factotLDLT_solvingFinalSync(FactorLDLTWorkerContext &ref_workerContext);
 307
 308     static int factotLDLT_scalingAndFactorizingFinal_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 309     static void factotLDLT_scalingAndFactorizingFinal(FactorLDLTWorkerContext &ref_workerContext, unsigned ownThreadIndex);
 310
 311     static int factotLDLT_completion_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 312
 313 private:
 314     struct FactorizationSolveL1StripeCellContext
 315     {
 316         template<unsigned int block_step, unsigned int b_rows>
 317         static void initializePrecalculatedZs(dReal (&Z)[block_step][b_rows])
 318         {
 319             Z[0][0] = 0;
 320             if (b_rows >= 2)
 321             {
 322                 Z[0][1] = 0;
 323             }
 324             Z[1][0] = 0;
 325             if (b_rows >= 2)
 326             {
 327                 Z[1][1] = 0;
 328             }
 329             dSASSERT(block_step == 2);
 330             dSASSERT(b_rows >= 1 && b_rows <= 2);
 331         }
 332
 333         template<unsigned int block_step, unsigned int b_rows>
 334         void loadPrecalculatedZs(dReal (&Z)[block_step][b_rows]) const
 335         {
 336             dSASSERT(block_step <= dARRAY_SIZE(m_c));
 337             dSASSERT(b_rows <= dARRAY_SIZE(m_c[0]));
 338
 339             Z[0][0] = m_c[0][0];
 340             if (b_rows >= 2)
 341             {
 342                 Z[0][1] = m_c[0][1];
 343             }
 344             Z[1][0] = m_c[1][0];
 345             if (b_rows >= 2)
 346             {
 347                 Z[1][1] = m_c[1][1];
 348             }
 349             dSASSERT(block_step == 2);
 350             dSASSERT(b_rows >= 1 && b_rows <= 2);
 351         }
 352
 353         template<unsigned int block_step, unsigned int b_rows>
 354         void storePrecalculatedZs(const dReal (&Z)[block_step][b_rows])
 355         {
 356             dSASSERT(block_step <= dARRAY_SIZE(m_c));
 357             dSASSERT(b_rows <= dARRAY_SIZE(m_c[0]));
 358
 359             m_c[0][0] = Z[0][0];
 360             if (b_rows >= 2)
 361             {
 362                 m_c[0][1] = Z[0][1];
 363             }
 364             m_c[1][0] = Z[1][0];
 365             if (b_rows >= 2)
 366             {
 367                 m_c[1][1] = Z[1][1];
 368             }
 369             dSASSERT(block_step == 2);
 370             dSASSERT(b_rows >= 1 && b_rows <= 2);
 371         }
 372
 373         dReal m_c[FSL1S_BLOCK_SIZE][FSL1S_REGULAR_B_ROWS];
 374         // dReal m_reserved[4];
 375     };
 376
 377     static FactorizationSolveL1StripeCellContext &buildBlockContextRef(FactorizationSolveL1StripeCellContext *cellContexts, unsigned blockIndex, CellContextInstance contextInstance)
 378     {
 379         return cellContexts[blockIndex * CCI__MAX + contextInstance];
 380     }
 381
 382     static FactorizationSolveL1StripeCellContext &buildResultContextRef(FactorizationSolveL1StripeCellContext *cellContexts, unsigned blockIndex, unsigned blockCount)
 383     {
 384         return cellContexts[blockCount * CCI__MAX + blockIndex];
 385     }
 386
 387 private:
 388     struct FactorizationFactorizeL1StripeThreadContext
 389     {
 390         template<unsigned int a_rows>
 391         void assignDataSum(const dReal (&sameZ)[a_rows], const dReal (&mixedZ)[dMACRO_MAX(a_rows - 1, 1)],
 392             const FactorizationFactorizeL1StripeThreadContext &partialSumContext)
 393         {
 394             m_sameZ[0] = sameZ[0] + partialSumContext.m_sameZ[0];
 395             if (a_rows >= 2)
 396             {
 397                 m_sameZ[1] = sameZ[1] + partialSumContext.m_sameZ[1];
 398                 m_mixedZ[0] = mixedZ[0] + partialSumContext.m_mixedZ[0];
 399             }
 400         }
 401
 402         template<unsigned int a_rows>
 403         void assignDataAlone(const dReal (&sameZ)[a_rows], const dReal (&mixedZ)[dMACRO_MAX(a_rows - 1, 1)])
 404         {
 405             m_sameZ[0] = sameZ[0];
 406             if (a_rows >= 2)
 407             {
 408                 m_sameZ[1] = sameZ[1];
 409                 m_mixedZ[0] = mixedZ[0];
 410             }
 411         }
 412
 413         template<unsigned int a_rows>
 414         void retrieveData(dReal (&out_sameZ)[a_rows], dReal (&out_mixedZ)[dMACRO_MAX(a_rows - 1, 1)]) const
 415         {
 416             out_sameZ[0] = m_sameZ[0];
 417             if (a_rows >= 2)
 418             {
 419                 out_sameZ[1] = m_sameZ[1];
 420                 out_mixedZ[0] = m_mixedZ[0];
 421             }
 422             dAASSERT(a_rows >= 1 && a_rows <= 2);
 423         }
 424
 425         dReal m_sameZ[FFL1S_REGULAR_A_ROWS];
 426         dReal m_mixedZ[dMACRO_MAX(FFL1S_REGULAR_A_ROWS - 1, 1)];
 427         dReal m_reserved[1]; // [5]; // for alignment
 428     };
 429
 430     struct FactorizationFactorizeL1StripeContext
 431     {
 432         void initialize(unsigned threadCount)
 433         {
 434             m_threadsRunning = threadCount;
 435             m_nextColumnIndex = 0;
 436             m_sumThreadIndex = 0;
 437         }
 438
 439         atomicord32 m_threadsRunning;
 440         atomicord32 m_nextColumnIndex;
 441         volatile atomicord32 m_sumThreadIndex;
 442         atomicord32 m_reserved[1]; // [13]; // for alignment
 443         FactorizationFactorizeL1StripeThreadContext m_threadContexts[1]; // =[threadCount]
 444     };
 445
 446 private:
 447     struct SolveL1StraightCellContext;
 448
 449     enum
 450     {
 451         SL1S_COOPERATIVE_BLOCK_COUNT_MINIMUM = 8,
 452
 453         SL1S_B_STRIDE   = 1,
 454         SL1S_BLOCK_SIZE = 4,
 455     };
 456
 457     static unsigned restrictSolvingL1StraightAllowedThreadCount(
 458         dxThreadingBase *threading, unsigned allowedThreadCount, unsigned rowCount);
 459     static void doEstimateCooperativeSolvingL1StraightResourceRequirementsValidated(
 460         dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
 461         unsigned allowedThreadCount, unsigned rowCount);
 462     static void doCooperativelySolveL1StraightValidated(
 463         dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
 464         const dReal *L, dReal *b, unsigned rowCount, unsigned rowSkip);
 465
 466     static unsigned deriveSolvingL1StraightBlockCount(unsigned rowCount, unsigned blockStep)
 467     {
 468         return (rowCount + (blockStep - 1)) / blockStep;
 469     }
 470
 471     struct SolvingL1StraightMemoryEstimates
 472     {
 473         void assignData(sizeint descriptorSizeRequired, sizeint contextSizeRequired)
 474         {
 475             m_descriptorSizeRequired = descriptorSizeRequired;
 476             m_contextSizeRequired = contextSizeRequired;
 477         }
 478
 479         sizeint  m_descriptorSizeRequired;
 480         sizeint  m_contextSizeRequired;
 481     };
 482
 483     static unsigned deriveSolvingL1StraightThreadCount(unsigned blockCount, unsigned allowedThreadCount)
 484     {
 485         dIASSERT(allowedThreadCount >= 1);
 486
 487         unsigned maximumCount = 1 + blockCount / SL1S_COOPERATIVE_BLOCK_COUNT_MINIMUM;
 488         return maximumCount >= allowedThreadCount ? allowedThreadCount : dMACRO_MAX(maximumCount, 1U);
 489     }
 490
 491     template<unsigned int block_step>
 492     static sizeint estimateCooperativelySolvingL1StraightMemoryRequirement(unsigned rowCount, SolvingL1StraightMemoryEstimates &ref_solvingMemoryEstimates);
 493
 494     static void *markCooperativelySolvingL1StraightMemoryStructuresOut(void *buffer,
 495         const SolvingL1StraightMemoryEstimates &solvingMemoryEstimates,
 496         cellindexint *&out_blockProgressDescriptors, SolveL1StraightCellContext *&out_cellContexts)
 497     {
 498         void *currentLocation = buffer;
 499
 500         out_blockProgressDescriptors = (cellindexint *)currentLocation; currentLocation = (uint8 *)currentLocation + solvingMemoryEstimates.m_descriptorSizeRequired;
 501         out_cellContexts = (SolveL1StraightCellContext *)currentLocation; currentLocation = (uint8 *)currentLocation + solvingMemoryEstimates.m_contextSizeRequired;
 502         return currentLocation;
 503     }
 504
 505     template<unsigned int block_step>
 506     static void initializeCooperativelySolveL1StraightMemoryStructures(unsigned rowCount,
 507         atomicord32 &out_blockCompletionProgress, cellindexint *blockProgressDescriptors, SolveL1StraightCellContext *cellContexts);
 508     template<unsigned int block_step, unsigned int b_stride>
 509     static void participateSolvingL1Straight(const dReal *L, dReal *B, unsigned rowCount, unsigned rowSkip,
 510         volatile atomicord32 &refBlockCompletionProgress/*=0*/, volatile cellindexint *blockProgressDescriptors/*=[blockCount]*/,
 511         SolveL1StraightCellContext *cellContexts/*=[CCI__MAX x blockCount] + [blockCount]*/, unsigned ownThreadIndex);
 512
 513 private:
 514     struct SolveL1StraightWorkerContext
 515     {
 516         void init(const dReal *L, dReal *b, unsigned rowCount, unsigned rowSkip,
 517             atomicord32 &ref_blockCompletionProgress, cellindexint *blockProgressDescriptors, SolveL1StraightCellContext *cellContexts)
 518         {
 519             m_L = L;
 520             m_b = b;
 521             m_rowCount = rowCount;
 522             m_rowSkip = rowSkip;
 523             m_ptrBlockCompletionProgress = &ref_blockCompletionProgress;
 524             m_blockProgressDescriptors = blockProgressDescriptors;
 525             m_cellContexts = cellContexts;
 526         }
 527
 528         const dReal     *m_L;
 529         dReal           *m_b;
 530         unsigned        m_rowCount;
 531         unsigned        m_rowSkip;
 532         atomicord32     *m_ptrBlockCompletionProgress;
 533         cellindexint    *m_blockProgressDescriptors;
 534         SolveL1StraightCellContext *m_cellContexts;
 535     };
 536
 537     static int solveL1Straight_worker_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 538     static void solveL1Straight_worker(SolveL1StraightWorkerContext &ref_context, unsigned ownThreadIndex);
 539
 540     static int solveL1Straight_completion_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 541
 542 private:
 543     struct SolveL1StraightCellContext
 544     {
 545         template<unsigned int block_step>
 546         static void initializePrecalculatedZs(dReal (&Z)[block_step])
 547         {
 548             std::fill(Z, Z + block_step, REAL(0.0));
 549         }
 550
 551         template<unsigned int block_step>
 552         void loadPrecalculatedZs(dReal (&Z)[block_step]) const
 553         {
 554             dSASSERT(block_step <= dARRAY_SIZE(m_c));
 555
 556             std::copy(m_c, m_c + block_step, Z);
 557         }
 558
 559         template<unsigned int block_step>
 560         void storePrecalculatedZs(const dReal (&Z)[block_step])
 561         {
 562             dSASSERT(block_step <= dARRAY_SIZE(m_c));
 563
 564             std::copy(Z, Z + block_step, m_c);
 565         }
 566
 567         dReal m_c[SL1S_BLOCK_SIZE];
 568     };
 569
 570
 571     static SolveL1StraightCellContext &buildBlockContextRef(SolveL1StraightCellContext *cellContexts, unsigned blockIndex, CellContextInstance contextInstance)
 572     {
 573         return cellContexts[blockIndex * CCI__MAX + contextInstance];
 574     }
 575
 576     static SolveL1StraightCellContext &buildResultContextRef(SolveL1StraightCellContext *cellContexts, unsigned blockIndex, unsigned blockCount)
 577     {
 578         return cellContexts[blockCount * CCI__MAX + blockIndex];
 579     }
 580
 581
 582 private:
 583     struct SolveL1TransposedCellContext;
 584
 585     enum
 586     {
 587         SL1T_COOPERATIVE_BLOCK_COUNT_MINIMUM = SL1S_COOPERATIVE_BLOCK_COUNT_MINIMUM,
 588
 589         SL1T_B_STRIDE   = SL1S_B_STRIDE,
 590         SL1T_BLOCK_SIZE = 4,
 591     };
 592
 593     static unsigned restrictSolvingL1TransposedAllowedThreadCount(
 594         dxThreadingBase *threading, unsigned allowedThreadCount, unsigned rowCount);
 595     static void doEstimateCooperativeSolvingL1TransposedResourceRequirementsValidated(
 596         dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
 597         unsigned allowedThreadCount, unsigned rowCount);
 598     static void doCooperativelySolveL1TransposedValidated(
 599         dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
 600         const dReal *L, dReal *b, unsigned rowCount, unsigned rowSkip);
 601
 602     static unsigned deriveSolvingL1TransposedBlockCount(unsigned rowCount, unsigned blockStep)
 603     {
 604         return (rowCount + (blockStep - 1)) / blockStep;
 605     }
 606
 607     struct SolvingL1TransposedMemoryEstimates
 608     {
 609         void assignData(sizeint descriptorSizeRequired, sizeint contextSizeRequired)
 610         {
 611             m_descriptorSizeRequired = descriptorSizeRequired;
 612             m_contextSizeRequired = contextSizeRequired;
 613         }
 614
 615         sizeint  m_descriptorSizeRequired;
 616         sizeint  m_contextSizeRequired;
 617     };
 618
 619     static unsigned deriveSolvingL1TransposedThreadCount(unsigned blockCount, unsigned allowedThreadCount)
 620     {
 621         dSASSERT(SL1T_COOPERATIVE_BLOCK_COUNT_MINIMUM + 0 == SL1S_COOPERATIVE_BLOCK_COUNT_MINIMUM);
 622
 623         return deriveSolvingL1StraightThreadCount(blockCount, allowedThreadCount);
 624     }
 625
 626     template<unsigned int block_step>
 627     static sizeint estimateCooperativelySolvingL1TransposedMemoryRequirement(unsigned rowCount, SolvingL1TransposedMemoryEstimates &ref_solvingMemoryEstimates);
 628
 629     static void *markCooperativelySolvingL1TransposedMemoryStructuresOut(void *buffer,
 630         const SolvingL1TransposedMemoryEstimates &solvingMemoryEstimates,
 631         cellindexint *&out_blockProgressDescriptors, SolveL1TransposedCellContext *&out_cellContexts)
 632     {
 633         void *currentLocation = buffer;
 634
 635         out_blockProgressDescriptors = (cellindexint *)currentLocation; currentLocation = (uint8 *)currentLocation + solvingMemoryEstimates.m_descriptorSizeRequired;
 636         out_cellContexts = (SolveL1TransposedCellContext *)currentLocation; currentLocation = (uint8 *)currentLocation + solvingMemoryEstimates.m_contextSizeRequired;
 637         return currentLocation;
 638     }
 639
 640     template<unsigned int block_step>
 641     static void *allocateCooperativelySolveL1TransposedMemoryStructures(sizeint &out_sizeAllocated, unsigned rowCount,
 642         cellindexint *&out_blockProgressDescriptors, SolveL1TransposedCellContext *&out_cellContexts);
 643     template<unsigned int block_step>
 644     static void initializeCooperativelySolveL1TransposedMemoryStructures(unsigned rowCount,
 645         atomicord32 &out_blockCompletionProgress, cellindexint *blockProgressDescriptors, SolveL1TransposedCellContext *cellContexts);
 646     template<unsigned int block_step, unsigned int b_stride>
 647     static void participateSolvingL1Transposed(const dReal *L, dReal *B, unsigned rowCount, unsigned rowSkip,
 648         volatile atomicord32 &refBlockCompletionProgress/*=0*/, volatile cellindexint *blockProgressDescriptors/*=[blockCount]*/,
 649         SolveL1TransposedCellContext *cellContexts/*=[CCI__MAX x blockCount] + [blockCount]*/, unsigned ownThreadIndex);
 650
 651 private:
 652     struct SolveL1TransposedWorkerContext
 653     {
 654         void init(const dReal *L, dReal *b, unsigned rowCount, unsigned rowSkip,
 655             atomicord32 &ref_blockCompletionProgress, cellindexint *blockProgressDescriptors, SolveL1TransposedCellContext *cellContexts)
 656         {
 657             m_L = L;
 658             m_b = b;
 659             m_rowCount = rowCount;
 660             m_rowSkip = rowSkip;
 661             m_ptrBlockCompletionProgress = &ref_blockCompletionProgress;
 662             m_blockProgressDescriptors = blockProgressDescriptors;
 663             m_cellContexts = cellContexts;
 664         }
 665
 666         const dReal     *m_L;
 667         dReal           *m_b;
 668         unsigned        m_rowCount;
 669         unsigned        m_rowSkip;
 670         atomicord32     *m_ptrBlockCompletionProgress;
 671         cellindexint    *m_blockProgressDescriptors;
 672         SolveL1TransposedCellContext *m_cellContexts;
 673     };
 674
 675     static int solveL1Transposed_worker_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 676     static void solveL1Transposed_worker(SolveL1TransposedWorkerContext &ref_context, unsigned ownThreadIndex);
 677
 678     static int solveL1Transposed_completion_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 679
 680 private:
 681     struct SolveL1TransposedCellContext
 682     {
 683         template<unsigned int block_step>
 684         static void initializePrecalculatedZs(dReal (&Z)[block_step])
 685         {
 686             std::fill(Z, Z + block_step, REAL(0.0));
 687         }
 688
 689         template<unsigned int block_step>
 690         void loadPrecalculatedZs(dReal (&Z)[block_step]) const
 691         {
 692             dSASSERT(block_step <= dARRAY_SIZE(m_c));
 693
 694             std::copy(m_c, m_c + block_step, Z);
 695         }
 696
 697         template<unsigned int block_step>
 698         void storePrecalculatedZs(const dReal (&Z)[block_step])
 699         {
 700             dSASSERT(block_step <= dARRAY_SIZE(m_c));
 701
 702             std::copy(Z, Z + block_step, m_c);
 703         }
 704
 705         dReal m_c[SL1T_BLOCK_SIZE];
 706     };
 707
 708     static SolveL1TransposedCellContext &buildBlockContextRef(SolveL1TransposedCellContext *cellContexts, unsigned blockIndex, CellContextInstance contextInstance)
 709     {
 710         return cellContexts[blockIndex * CCI__MAX + contextInstance];
 711     }
 712
 713     static SolveL1TransposedCellContext &buildResultContextRef(SolveL1TransposedCellContext *cellContexts, unsigned blockIndex, unsigned blockCount)
 714     {
 715         return cellContexts[blockCount * CCI__MAX + blockIndex];
 716     }
 717
 718 private:
 719     enum
 720     {
 721         SV_A_STRIDE = 1,
 722         SV_D_STRIDE = 1,
 723
 724         SV_BLOCK_SIZE = 128,
 725         SV_COOPERATIVE_BLOCK_COUNT_MINIMUM = 3,
 726     };
 727
 728     static unsigned restrictScalingVectorAllowedThreadCount(
 729         dxThreadingBase *threading, unsigned allowedThreadCount, unsigned elementCount);
 730     static void doEstimateCooperativeScalingVectorResourceRequirementsValidated(
 731         dxResourceRequirementDescriptor *summaryRequirementsDescriptor,
 732         unsigned allowedThreadCount, unsigned elementCount);
 733     static void doCooperativelyScaleVectorValidated(dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount,
 734         dReal *vectorData, const dReal *scaleData, unsigned elementCount);
 735
 736     static unsigned deriveScalingVectorBlockCount(unsigned elementCount, unsigned blockStep)
 737     {
 738         return (elementCount + (blockStep - 1)) / blockStep;
 739     }
 740
 741     static unsigned deriveScalingVectorThreadCount(unsigned lastBlockIndex, unsigned allowedThreadCount)
 742     {
 743         dIASSERT(allowedThreadCount >= 1);
 744
 745         unsigned maximumCount = lastBlockIndex;
 746         return maximumCount >= allowedThreadCount ? allowedThreadCount : dMACRO_MAX(maximumCount, 1U);
 747     }
 748
 749     static void initializeCooperativelyScaleVectorMemoryStructures(atomicord32 &out_blockCompletionProgress)
 750     {
 751         out_blockCompletionProgress = 0;
 752     }
 753     template<unsigned int block_step, unsigned int a_stride, unsigned int d_stride>
 754     static void participateScalingVector(dReal *ptrAStart, const dReal *ptrDStart, const unsigned elementCount,
 755         volatile atomicord32 &refBlockCompletionProgress/*=0*/);
 756
 757 private:
 758     struct ScaleVectorWorkerContext
 759     {
 760         void init(dReal *vectorData, const dReal *scaleData, unsigned elementCount,
 761             atomicord32 &ref_blockCompletionProgress)
 762         {
 763             m_vectorData = vectorData;
 764             m_scaleData = scaleData;
 765             m_elementCount = elementCount;
 766             m_ptrBlockCompletionProgress = &ref_blockCompletionProgress;
 767         }
 768
 769         dReal           *m_vectorData;
 770         const dReal     *m_scaleData;
 771         unsigned        m_elementCount;
 772         atomicord32     *m_ptrBlockCompletionProgress;
 773     };
 774
 775     static int scaleVector_worker_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 776     static void scaleVector_worker(ScaleVectorWorkerContext &ref_context);
 777
 778     static int scaleVector_completion_callback(void *callContext, dcallindex_t callInstanceIndex, dCallReleaseeID callThisReleasee);
 779
 780
 781 private:
 782     enum SolvingLDLTStage
 783     {
 784         SLDLTS__MIN,
 785
 786         SLDLTS_SOLVING_STRAIGHT = SLDLTS__MIN,
 787         SLDLTS_SCALING_VECTOR,
 788         SLDLTS_SOLVING_TRANSPOSED,
 789
 790         SLDLTS__MAX,
 791     };
 792
 793     enum
 794     {
 795         SLDLT_B_STRIDE          = SL1S_B_STRIDE,
 796         SLDLT_D_STRIDE          = FLDLT_D_STRIDE,
 797     };
 798
 799     static unsigned restrictSolvingLDLTAllowedThreadCount(
 800         dxThreadingBase *threading, unsigned allowedThreadCount, unsigned rowCount, unsigned &out_stageBlockCountSifficiencyMask);
 801
 802     static void doCooperativelySolveLDLTValidated(
 803         dxRequiredResourceContainer *resourceContainer, unsigned allowedThreadCount, unsigned stageBlockCountSifficiencyMask,
 804         const dReal *L, const dReal *d, dReal *b, unsigned rowCount, unsigned rowSkip);
 805 };
 806
 807
 808 #endif
 809