/[escript]/branches/doubleplusgood/paso/src/BlockOps.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/BlockOps.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4261 - (show annotations)
Wed Feb 27 06:09:33 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 1781 byte(s)
Initial all c++ build.
But ... there are now reinterpret_cast<>'s
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 #include "BlockOps.h"
18 #include "Paso.h"
19
20 /************************************************************************************
21
22 PASO block operations
23
24 ************************************************************************************/
25
26 void Paso_BlockOps_solveAll(dim_t n_block,dim_t n,double* D,index_t* pivot,double* x) {
27 dim_t i;
28 int failed=0;
29 const dim_t block_size=n_block*n_block;
30 (void)block_size; /* silence warning from var being unused by macros */
31
32 if (n_block==1) {
33 #pragma omp parallel for private(i) schedule(static)
34 for (i=0;i<n;++i) x[i]*=D[i];
35 } else if (n_block==2) {
36 #pragma omp parallel for private(i) schedule(static)
37 for (i=0;i<n;++i) Paso_BlockOps_MViP_2(&D[4*i], &x[2*i]);
38
39 } else if (n_block==3) {
40 #pragma omp parallel for private(i) schedule(static)
41 for (i=0;i<n;++i) Paso_BlockOps_MViP_3(&D[9*i], &x[3*i]);
42 } else {
43
44 #pragma omp parallel for private(i) schedule(static)
45 for (i=0;i<n;++i) {
46 Paso_BlockOps_solve_N(n_block, &x[n_block*i], &D[block_size*i], &pivot[n_block*i], &failed);
47 }
48 }
49 if (failed > 0) {
50 Esys_setError(ZERO_DIVISION_ERROR, "Paso_BlockOps_solveAll: solution failed.");
51 }
52 return;
53 }

  ViewVC Help
Powered by ViewVC 1.1.26