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"><../src/mat/impls/baij/seq/baij.h></A>
<a name="line2"> 2: </a>#include <A href="../../../../include/petsc/private/kernels/blockinvert.h.html"><petsc/private/kernels/blockinvert.h></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->data;
<a name="line7"> 7: </a> const <a href="../../../../../manualpages/Sys/PetscInt.html">PetscInt</a> *diag = a->diag, n = a->mbs, *vi, *ai = a->i, *aj = a->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->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, &b));
<a name="line15"> 15: </a> <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecGetArray.html">VecGetArray</a>(xx, &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 < 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 >= 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, &b));
<a name="line87"> 87: </a> <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecRestoreArray.html">VecRestoreArray</a>(xx, &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->nz) - 5.0 * A->cmap->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->data;
<a name="line95"> 95: </a> const <a href="../../../../../manualpages/Sys/PetscInt.html">PetscInt</a> n = a->mbs, *vi, *ai = a->i, *aj = a->j, *adiag = a->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->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, &b));
<a name="line103">103: </a> <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecGetArray.html">VecGetArray</a>(xx, &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 < 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 < 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 >= 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 < 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, &b));
<a name="line176">176: </a> <a href="../../../../../manualpages/Sys/PetscCall.html">PetscCall</a>(<a href="../../../../../manualpages/Vec/VecRestoreArray.html">VecRestoreArray</a>(xx, &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->nz) - 5.0 * A->cmap->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>
|