File: baijsolvnat5.c.html

package info (click to toggle)
petsc 3.22.5%2Bdfsg1-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 516,740 kB
  • sloc: ansic: 814,333; cpp: 50,948; python: 37,416; f90: 17,187; javascript: 3,493; makefile: 3,198; sh: 1,502; xml: 619; objc: 445; java: 13; csh: 1
file content (195 lines) | stat: -rw-r--r-- 14,472 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
<center><a href="https://gitlab.com/petsc/petsc/-/blob/9fe822445bcdd45fb93170ff68fa7403d3f52f09/src/mat/impls/baij/seq/baijsolvnat5.c">Actual source code: baijsolvnat5.c</a></center><br>

<html>
<head>
<title></title>
<meta name="generator" content="c2html 0.9.6">
<meta name="date" content="2025-03-28T21:11:00+00:00">
</head>

<body bgcolor="#FFFFFF">
<pre width=80>
<a name="line1">  1: </a>#include <A href="../../../../src/mat/impls/baij/seq/baij.h.html">&lt;../src/mat/impls/baij/seq/baij.h&gt;</A>
<a name="line2">  2: </a>#include <A href="../../../../include/petsc/private/kernels/blockinvert.h.html">&lt;petsc/private/kernels/blockinvert.h&gt;</A>

<a name="line4">  4: </a><strong><font color="#4169E1"><a name="MatSolve_SeqBAIJ_5_NaturalOrdering_inplace"></a><a href="../../../../../manualpages/Sys/PetscErrorCode.html">PetscErrorCode</a> MatSolve_SeqBAIJ_5_NaturalOrdering_inplace(<a href="../../../../../manualpages/Mat/Mat.html">Mat</a> A, <a href="../../../../../manualpages/Vec/Vec.html">Vec</a> bb, <a href="../../../../../manualpages/Vec/Vec.html">Vec</a> xx)</font></strong>
<a name="line5">  5: </a>{
<a name="line6">  6: </a>  Mat_SeqBAIJ       *a    = (Mat_SeqBAIJ *)A-&gt;data;
<a name="line7">  7: </a>  const <a href="../../../../../manualpages/Sys/PetscInt.html">PetscInt</a>    *diag = a-&gt;diag, n = a-&gt;mbs, *vi, *ai = a-&gt;i, *aj = a-&gt;j;
<a name="line8">  8: </a>  <a href="../../../../../manualpages/Sys/PetscInt.html">PetscInt</a>           i, nz, idx, idt, jdx;
<a name="line9">  9: </a>  const MatScalar   *aa = a-&gt;a, *v;
<a name="line10"> 10: </a>  <a href="../../../../../manualpages/Sys/PetscScalar.html">PetscScalar</a>       *x, s1, s2, s3, s4, s5, x1, x2, x3, x4, x5;
<a name="line11"> 11: </a>  const <a href="../../../../../manualpages/Sys/PetscScalar.html">PetscScalar</a> *b;

<a name="line13"> 13: </a>  <a href="../../../../../manualpages/Sys/PetscFunctionBegin.html">PetscFunctionBegin</a>;
<a name="line14"> 14: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecGetArrayRead.html">VecGetArrayRead</a>(bb, &amp;b));
<a name="line15"> 15: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecGetArray.html">VecGetArray</a>(xx, &amp;x));
<a name="line16"> 16: </a>  <font color="#B22222">/* forward solve the lower triangular */</font>
<a name="line17"> 17: </a>  idx  = 0;
<a name="line18"> 18: </a>  x[0] = b[idx];
<a name="line19"> 19: </a>  x[1] = b[1 + idx];
<a name="line20"> 20: </a>  x[2] = b[2 + idx];
<a name="line21"> 21: </a>  x[3] = b[3 + idx];
<a name="line22"> 22: </a>  x[4] = b[4 + idx];
<a name="line23"> 23: </a>  <font color="#4169E1">for</font> (i = 1; i &lt; n; i++) {
<a name="line24"> 24: </a>    v   = aa + 25 * ai[i];
<a name="line25"> 25: </a>    vi  = aj + ai[i];
<a name="line26"> 26: </a>    nz  = diag[i] - ai[i];
<a name="line27"> 27: </a>    idx = 5 * i;
<a name="line28"> 28: </a>    s1  = b[idx];
<a name="line29"> 29: </a>    s2  = b[1 + idx];
<a name="line30"> 30: </a>    s3  = b[2 + idx];
<a name="line31"> 31: </a>    s4  = b[3 + idx];
<a name="line32"> 32: </a>    s5  = b[4 + idx];
<a name="line33"> 33: </a>    <font color="#4169E1">while</font> (nz--) {
<a name="line34"> 34: </a>      jdx = 5 * (*vi++);
<a name="line35"> 35: </a>      x1  = x[jdx];
<a name="line36"> 36: </a>      x2  = x[1 + jdx];
<a name="line37"> 37: </a>      x3  = x[2 + jdx];
<a name="line38"> 38: </a>      x4  = x[3 + jdx];
<a name="line39"> 39: </a>      x5  = x[4 + jdx];
<a name="line40"> 40: </a>      s1 -= v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
<a name="line41"> 41: </a>      s2 -= v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
<a name="line42"> 42: </a>      s3 -= v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
<a name="line43"> 43: </a>      s4 -= v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
<a name="line44"> 44: </a>      s5 -= v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
<a name="line45"> 45: </a>      v += 25;
<a name="line46"> 46: </a>    }
<a name="line47"> 47: </a>    x[idx]     = s1;
<a name="line48"> 48: </a>    x[1 + idx] = s2;
<a name="line49"> 49: </a>    x[2 + idx] = s3;
<a name="line50"> 50: </a>    x[3 + idx] = s4;
<a name="line51"> 51: </a>    x[4 + idx] = s5;
<a name="line52"> 52: </a>  }
<a name="line53"> 53: </a>  <font color="#B22222">/* backward solve the upper triangular */</font>
<a name="line54"> 54: </a>  <font color="#4169E1">for</font> (i = n - 1; i &gt;= 0; i--) {
<a name="line55"> 55: </a>    v   = aa + 25 * diag[i] + 25;
<a name="line56"> 56: </a>    vi  = aj + diag[i] + 1;
<a name="line57"> 57: </a>    nz  = ai[i + 1] - diag[i] - 1;
<a name="line58"> 58: </a>    idt = 5 * i;
<a name="line59"> 59: </a>    s1  = x[idt];
<a name="line60"> 60: </a>    s2  = x[1 + idt];
<a name="line61"> 61: </a>    s3  = x[2 + idt];
<a name="line62"> 62: </a>    s4  = x[3 + idt];
<a name="line63"> 63: </a>    s5  = x[4 + idt];
<a name="line64"> 64: </a>    <font color="#4169E1">while</font> (nz--) {
<a name="line65"> 65: </a>      idx = 5 * (*vi++);
<a name="line66"> 66: </a>      x1  = x[idx];
<a name="line67"> 67: </a>      x2  = x[1 + idx];
<a name="line68"> 68: </a>      x3  = x[2 + idx];
<a name="line69"> 69: </a>      x4  = x[3 + idx];
<a name="line70"> 70: </a>      x5  = x[4 + idx];
<a name="line71"> 71: </a>      s1 -= v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
<a name="line72"> 72: </a>      s2 -= v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
<a name="line73"> 73: </a>      s3 -= v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
<a name="line74"> 74: </a>      s4 -= v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
<a name="line75"> 75: </a>      s5 -= v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
<a name="line76"> 76: </a>      v += 25;
<a name="line77"> 77: </a>    }
<a name="line78"> 78: </a>    v          = aa + 25 * diag[i];
<a name="line79"> 79: </a>    x[idt]     = v[0] * s1 + v[5] * s2 + v[10] * s3 + v[15] * s4 + v[20] * s5;
<a name="line80"> 80: </a>    x[1 + idt] = v[1] * s1 + v[6] * s2 + v[11] * s3 + v[16] * s4 + v[21] * s5;
<a name="line81"> 81: </a>    x[2 + idt] = v[2] * s1 + v[7] * s2 + v[12] * s3 + v[17] * s4 + v[22] * s5;
<a name="line82"> 82: </a>    x[3 + idt] = v[3] * s1 + v[8] * s2 + v[13] * s3 + v[18] * s4 + v[23] * s5;
<a name="line83"> 83: </a>    x[4 + idt] = v[4] * s1 + v[9] * s2 + v[14] * s3 + v[19] * s4 + v[24] * s5;
<a name="line84"> 84: </a>  }

<a name="line86"> 86: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecRestoreArrayRead.html">VecRestoreArrayRead</a>(bb, &amp;b));
<a name="line87"> 87: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecRestoreArray.html">VecRestoreArray</a>(xx, &amp;x));
<a name="line88"> 88: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Profiling/PetscLogFlops.html">PetscLogFlops</a>(2.0 * 25 * (a-&gt;nz) - 5.0 * A-&gt;cmap-&gt;n));
<a name="line89"> 89: </a>  <a href="../../../../../manualpages/Sys/PetscFunctionReturn.html">PetscFunctionReturn</a>(<a href="../../../../../manualpages/Sys/PetscErrorCode.html">PETSC_SUCCESS</a>);
<a name="line90"> 90: </a>}

<a name="line92"> 92: </a><strong><font color="#4169E1"><a name="MatSolve_SeqBAIJ_5_NaturalOrdering"></a><a href="../../../../../manualpages/Sys/PetscErrorCode.html">PetscErrorCode</a> MatSolve_SeqBAIJ_5_NaturalOrdering(<a href="../../../../../manualpages/Mat/Mat.html">Mat</a> A, <a href="../../../../../manualpages/Vec/Vec.html">Vec</a> bb, <a href="../../../../../manualpages/Vec/Vec.html">Vec</a> xx)</font></strong>
<a name="line93"> 93: </a>{
<a name="line94"> 94: </a>  Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A-&gt;data;
<a name="line95"> 95: </a>  const <a href="../../../../../manualpages/Sys/PetscInt.html">PetscInt</a>     n = a-&gt;mbs, *vi, *ai = a-&gt;i, *aj = a-&gt;j, *adiag = a-&gt;diag;
<a name="line96"> 96: </a>  <a href="../../../../../manualpages/Sys/PetscInt.html">PetscInt</a>           i, k, nz, idx, idt, jdx;
<a name="line97"> 97: </a>  const MatScalar   *aa = a-&gt;a, *v;
<a name="line98"> 98: </a>  <a href="../../../../../manualpages/Sys/PetscScalar.html">PetscScalar</a>       *x, s1, s2, s3, s4, s5, x1, x2, x3, x4, x5;
<a name="line99"> 99: </a>  const <a href="../../../../../manualpages/Sys/PetscScalar.html">PetscScalar</a> *b;

<a name="line101">101: </a>  <a href="../../../../../manualpages/Sys/PetscFunctionBegin.html">PetscFunctionBegin</a>;
<a name="line102">102: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecGetArrayRead.html">VecGetArrayRead</a>(bb, &amp;b));
<a name="line103">103: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecGetArray.html">VecGetArray</a>(xx, &amp;x));
<a name="line104">104: </a>  <font color="#B22222">/* forward solve the lower triangular */</font>
<a name="line105">105: </a>  idx  = 0;
<a name="line106">106: </a>  x[0] = b[idx];
<a name="line107">107: </a>  x[1] = b[1 + idx];
<a name="line108">108: </a>  x[2] = b[2 + idx];
<a name="line109">109: </a>  x[3] = b[3 + idx];
<a name="line110">110: </a>  x[4] = b[4 + idx];
<a name="line111">111: </a>  <font color="#4169E1">for</font> (i = 1; i &lt; n; i++) {
<a name="line112">112: </a>    v   = aa + 25 * ai[i];
<a name="line113">113: </a>    vi  = aj + ai[i];
<a name="line114">114: </a>    nz  = ai[i + 1] - ai[i];
<a name="line115">115: </a>    idx = 5 * i;
<a name="line116">116: </a>    s1  = b[idx];
<a name="line117">117: </a>    s2  = b[1 + idx];
<a name="line118">118: </a>    s3  = b[2 + idx];
<a name="line119">119: </a>    s4  = b[3 + idx];
<a name="line120">120: </a>    s5  = b[4 + idx];
<a name="line121">121: </a>    <font color="#4169E1">for</font> (k = 0; k &lt; nz; k++) {
<a name="line122">122: </a>      jdx = 5 * vi[k];
<a name="line123">123: </a>      x1  = x[jdx];
<a name="line124">124: </a>      x2  = x[1 + jdx];
<a name="line125">125: </a>      x3  = x[2 + jdx];
<a name="line126">126: </a>      x4  = x[3 + jdx];
<a name="line127">127: </a>      x5  = x[4 + jdx];
<a name="line128">128: </a>      s1 -= v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
<a name="line129">129: </a>      s2 -= v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
<a name="line130">130: </a>      s3 -= v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
<a name="line131">131: </a>      s4 -= v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
<a name="line132">132: </a>      s5 -= v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
<a name="line133">133: </a>      v += 25;
<a name="line134">134: </a>    }
<a name="line135">135: </a>    x[idx]     = s1;
<a name="line136">136: </a>    x[1 + idx] = s2;
<a name="line137">137: </a>    x[2 + idx] = s3;
<a name="line138">138: </a>    x[3 + idx] = s4;
<a name="line139">139: </a>    x[4 + idx] = s5;
<a name="line140">140: </a>  }

<a name="line142">142: </a>  <font color="#B22222">/* backward solve the upper triangular */</font>
<a name="line143">143: </a>  <font color="#4169E1">for</font> (i = n - 1; i &gt;= 0; i--) {
<a name="line144">144: </a>    v   = aa + 25 * (adiag[i + 1] + 1);
<a name="line145">145: </a>    vi  = aj + adiag[i + 1] + 1;
<a name="line146">146: </a>    nz  = adiag[i] - adiag[i + 1] - 1;
<a name="line147">147: </a>    idt = 5 * i;
<a name="line148">148: </a>    s1  = x[idt];
<a name="line149">149: </a>    s2  = x[1 + idt];
<a name="line150">150: </a>    s3  = x[2 + idt];
<a name="line151">151: </a>    s4  = x[3 + idt];
<a name="line152">152: </a>    s5  = x[4 + idt];
<a name="line153">153: </a>    <font color="#4169E1">for</font> (k = 0; k &lt; nz; k++) {
<a name="line154">154: </a>      idx = 5 * vi[k];
<a name="line155">155: </a>      x1  = x[idx];
<a name="line156">156: </a>      x2  = x[1 + idx];
<a name="line157">157: </a>      x3  = x[2 + idx];
<a name="line158">158: </a>      x4  = x[3 + idx];
<a name="line159">159: </a>      x5  = x[4 + idx];
<a name="line160">160: </a>      s1 -= v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
<a name="line161">161: </a>      s2 -= v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
<a name="line162">162: </a>      s3 -= v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
<a name="line163">163: </a>      s4 -= v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
<a name="line164">164: </a>      s5 -= v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
<a name="line165">165: </a>      v += 25;
<a name="line166">166: </a>    }
<a name="line167">167: </a>    <font color="#B22222">/* x = inv_diagonal*x */</font>
<a name="line168">168: </a>    x[idt]     = v[0] * s1 + v[5] * s2 + v[10] * s3 + v[15] * s4 + v[20] * s5;
<a name="line169">169: </a>    x[1 + idt] = v[1] * s1 + v[6] * s2 + v[11] * s3 + v[16] * s4 + v[21] * s5;
<a name="line170">170: </a>    x[2 + idt] = v[2] * s1 + v[7] * s2 + v[12] * s3 + v[17] * s4 + v[22] * s5;
<a name="line171">171: </a>    x[3 + idt] = v[3] * s1 + v[8] * s2 + v[13] * s3 + v[18] * s4 + v[23] * s5;
<a name="line172">172: </a>    x[4 + idt] = v[4] * s1 + v[9] * s2 + v[14] * s3 + v[19] * s4 + v[24] * s5;
<a name="line173">173: </a>  }

<a name="line175">175: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecRestoreArrayRead.html">VecRestoreArrayRead</a>(bb, &amp;b));
<a name="line176">176: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecRestoreArray.html">VecRestoreArray</a>(xx, &amp;x));
<a name="line177">177: </a>  <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Profiling/PetscLogFlops.html">PetscLogFlops</a>(2.0 * 25 * (a-&gt;nz) - 5.0 * A-&gt;cmap-&gt;n));
<a name="line178">178: </a>  <a href="../../../../../manualpages/Sys/PetscFunctionReturn.html">PetscFunctionReturn</a>(<a href="../../../../../manualpages/Sys/PetscErrorCode.html">PETSC_SUCCESS</a>);
<a name="line179">179: </a>}
</pre>
</body>

</html>