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

Contents of /branches/doubleplusgood/paso/src/Pattern_multiply.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: 3167 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
18 /************************************************************************************/
19
20 /* Paso: Pattern */
21
22 /************************************************************************************/
23
24 /* Author: Lutz Gross, l.gross@uq.edu.au */
25
26 /************************************************************************************/
27
28 #include "Paso.h"
29 #include "Pattern.h"
30
31 /************************************************************************************/
32
33 /* computes the pattern coming from matrix-matrix multiplication
34 *
35 **/
36
37 Paso_Pattern* Paso_Pattern_multiply(int type, Paso_Pattern* A, Paso_Pattern* B) {
38 Paso_Pattern*out=NULL;
39 index_t iptrA,iptrB;
40 dim_t i,j,k;
41 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(A->numOutput);
42
43 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
44 for(i = 0; i < A->numOutput; i++) {
45 for(iptrA = A->ptr[i]; iptrA < A->ptr[i+1]; ++iptrA) {
46 j = A->index[iptrA];
47 for(iptrB = B->ptr[j]; iptrB < B->ptr[j+1]; ++iptrB) {
48 k = B->index[iptrB];
49 Paso_IndexListArray_insertIndex(index_list,i,k);
50 }
51 }
52 }
53
54 out=Paso_Pattern_fromIndexListArray(0, index_list,0,B->numInput,0);
55
56 /* clean up */
57 Paso_IndexListArray_free(index_list);
58
59 return out;
60 }
61
62
63
64 /*
65 * Computes the pattern of C = A binary operation B for CSR matrices A,B
66 *
67 * Note: we do not check whether A_ij(op)B_ij=0
68 *
69 */
70 Paso_Pattern* Paso_Pattern_binop(int type, Paso_Pattern* A, Paso_Pattern* B) {
71 Paso_Pattern *out=NULL;
72 index_t iptrA,iptrB;
73 dim_t i,j,k;
74
75 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(A->numOutput);
76
77 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
78 for(i = 0; i < B->numOutput; i++){
79 iptrA = A->ptr[i],
80 iptrB = B->ptr[i];
81
82 while (iptrA < A->ptr[i+1] && iptrB < B->ptr[i+1]) {
83 j = A->index[iptrA];
84 k = B->index[iptrB];
85 if (j<k) {
86 Paso_IndexListArray_insertIndex(index_list,i,j);
87 iptrA++;
88 } else if (j>k) {
89 Paso_IndexListArray_insertIndex(index_list,i,k);
90 iptrB++;
91 } else if (j==k) {
92 Paso_IndexListArray_insertIndex(index_list,i,j);
93 iptrB++;
94 iptrA++;
95 }
96 }
97 while(iptrA < A->ptr[i+1]) {
98 j = A->index[iptrA];
99 Paso_IndexListArray_insertIndex(index_list,i,j);
100 iptrA++;
101 }
102 while(iptrB < B->ptr[i+1]) {
103 k = B->index[iptrB];
104 Paso_IndexListArray_insertIndex(index_list,i,k);
105 iptrB++;
106 }
107 }
108
109 out=Paso_Pattern_fromIndexListArray(0, index_list,0,A->numInput,0);
110
111
112 /* clean up */
113 Paso_IndexListArray_free(index_list);
114
115 return out;
116 }

  ViewVC Help
Powered by ViewVC 1.1.26