From: Martin Kronbichler Date: Mon, 2 Oct 2023 06:49:20 +0000 (+0200) Subject: Restructure FullMatrix::invert to (hopefully) lower compilation times X-Git-Tag: relicensing~445^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16076%2Fhead;p=dealii.git Restructure FullMatrix::invert to (hopefully) lower compilation times --- diff --git a/include/deal.II/lac/full_matrix.templates.h b/include/deal.II/lac/full_matrix.templates.h index 515ed92362..4f8163a4e9 100644 --- a/include/deal.II/lac/full_matrix.templates.h +++ b/include/deal.II/lac/full_matrix.templates.h @@ -1342,182 +1342,141 @@ FullMatrix::invert(const FullMatrix &M) Assert(this->n_rows() == M.n_rows(), ExcDimensionMismatch(this->n_rows(), M.n_rows())); - if (PointerComparison::equal(&M, this)) + switch (this->n_cols()) { - // avoid overwriting source - // by destination matrix: - const FullMatrix M2 = *this; - invert(M2); - } - else - switch (this->n_cols()) - { - case 1: - (*this)(0, 0) = number2(1.0) / M(0, 0); + case 1: + (*this)(0, 0) = number2(1.0) / M(0, 0); + return; + case 2: + { + const number2 M00 = M(0, 0); + const number2 M01 = M(0, 1); + const number2 M10 = M(1, 0); + const number2 M11 = M(1, 1); + const number2 t4 = number2(1.0) / (M00 * M11 - M01 * M10); + (*this)(0, 0) = M11 * t4; + (*this)(0, 1) = -M01 * t4; + (*this)(1, 0) = -M10 * t4; + (*this)(1, 1) = M00 * t4; + return; + }; + + case 3: + { + const number2 M00 = M(0, 0); + const number2 M01 = M(0, 1); + const number2 M02 = M(0, 2); + const number2 M10 = M(1, 0); + const number2 M11 = M(1, 1); + const number2 M12 = M(1, 2); + const number2 M20 = M(2, 0); + const number2 M21 = M(2, 1); + const number2 M22 = M(2, 2); + const number2 t00 = M11 * M22 - M12 * M21; + const number2 t10 = M12 * M20 - M10 * M22; + const number2 t20 = M10 * M21 - M11 * M20; + const number2 inv_det = + number2(1.0) / (M00 * t00 + M01 * t10 + M02 * t20); + (*this)(0, 0) = t00 * inv_det; + (*this)(0, 1) = (M02 * M21 - M01 * M22) * inv_det; + (*this)(0, 2) = (M01 * M12 - M02 * M11) * inv_det; + (*this)(1, 0) = t10 * inv_det; + (*this)(1, 1) = (M00 * M22 - M02 * M20) * inv_det; + (*this)(1, 2) = (M02 * M10 - M00 * M12) * inv_det; + (*this)(2, 0) = t20 * inv_det; + (*this)(2, 1) = (M01 * M20 - M00 * M21) * inv_det; + (*this)(2, 2) = (M00 * M11 - M01 * M10) * inv_det; return; - case 2: - // this is Maple output, - // thus a bit unstructured - { - const number2 t4 = - number2(1.0) / (M(0, 0) * M(1, 1) - M(0, 1) * M(1, 0)); - (*this)(0, 0) = M(1, 1) * t4; - (*this)(0, 1) = -M(0, 1) * t4; - (*this)(1, 0) = -M(1, 0) * t4; - (*this)(1, 1) = M(0, 0) * t4; - return; - }; - - case 3: - { - const number2 t4 = M(0, 0) * M(1, 1), t6 = M(0, 0) * M(1, 2), - t8 = M(0, 1) * M(1, 0), t00 = M(0, 2) * M(1, 0), - t01 = M(0, 1) * M(2, 0), t04 = M(0, 2) * M(2, 0), - t07 = number2(1.0) / - (t4 * M(2, 2) - t6 * M(2, 1) - t8 * M(2, 2) + - t00 * M(2, 1) + t01 * M(1, 2) - t04 * M(1, 1)); - (*this)(0, 0) = (M(1, 1) * M(2, 2) - M(1, 2) * M(2, 1)) * t07; - (*this)(0, 1) = -(M(0, 1) * M(2, 2) - M(0, 2) * M(2, 1)) * t07; - (*this)(0, 2) = -(-M(0, 1) * M(1, 2) + M(0, 2) * M(1, 1)) * t07; - (*this)(1, 0) = -(M(1, 0) * M(2, 2) - M(1, 2) * M(2, 0)) * t07; - (*this)(1, 1) = (M(0, 0) * M(2, 2) - t04) * t07; - (*this)(1, 2) = -(t6 - t00) * t07; - (*this)(2, 0) = -(-M(1, 0) * M(2, 1) + M(1, 1) * M(2, 0)) * t07; - (*this)(2, 1) = -(M(0, 0) * M(2, 1) - t01) * t07; - (*this)(2, 2) = (t4 - t8) * t07; - return; - }; - - case 4: - { - // with (linalg); - // a:=matrix(4,4); - // evalm(a); - // ai:=inverse(a); - // readlib(C); - // C(ai,optimized,filename=x4); - - const number2 t14 = M(0, 0) * M(1, 1); - const number2 t15 = M(2, 2) * M(3, 3); - const number2 t17 = M(2, 3) * M(3, 2); - const number2 t19 = M(0, 0) * M(2, 1); - const number2 t20 = M(1, 2) * M(3, 3); - const number2 t22 = M(1, 3) * M(3, 2); - const number2 t24 = M(0, 0) * M(3, 1); - const number2 t25 = M(1, 2) * M(2, 3); - const number2 t27 = M(1, 3) * M(2, 2); - const number2 t29 = M(1, 0) * M(0, 1); - const number2 t32 = M(1, 0) * M(2, 1); - const number2 t33 = M(0, 2) * M(3, 3); - const number2 t35 = M(0, 3) * M(3, 2); - const number2 t37 = M(1, 0) * M(3, 1); - const number2 t38 = M(0, 2) * M(2, 3); - const number2 t40 = M(0, 3) * M(2, 2); - const number2 t42 = t14 * t15 - t14 * t17 - t19 * t20 + t19 * t22 + - t24 * t25 - t24 * t27 - t29 * t15 + t29 * t17 + - t32 * t33 - t32 * t35 - t37 * t38 + t37 * t40; - const number2 t43 = M(2, 0) * M(0, 1); - const number2 t46 = M(2, 0) * M(1, 1); - const number2 t49 = M(2, 0) * M(3, 1); - const number2 t50 = M(0, 2) * M(1, 3); - const number2 t52 = M(0, 3) * M(1, 2); - const number2 t54 = M(3, 0) * M(0, 1); - const number2 t57 = M(3, 0) * M(1, 1); - const number2 t60 = M(3, 0) * M(2, 1); - const number2 t63 = t43 * t20 - t43 * t22 - t46 * t33 + t46 * t35 + - t49 * t50 - t49 * t52 - t54 * t25 + t54 * t27 + - t57 * t38 - t57 * t40 - t60 * t50 + t60 * t52; - const number2 t65 = number2(1.) / (t42 + t63); - const number2 t71 = M(0, 2) * M(2, 1); - const number2 t73 = M(0, 3) * M(2, 1); - const number2 t75 = M(0, 2) * M(3, 1); - const number2 t77 = M(0, 3) * M(3, 1); - const number2 t81 = M(0, 1) * M(1, 2); - const number2 t83 = M(0, 1) * M(1, 3); - const number2 t85 = M(0, 2) * M(1, 1); - const number2 t87 = M(0, 3) * M(1, 1); - const number2 t101 = M(1, 0) * M(2, 2); - const number2 t103 = M(1, 0) * M(2, 3); - const number2 t105 = M(2, 0) * M(1, 2); - const number2 t107 = M(2, 0) * M(1, 3); - const number2 t109 = M(3, 0) * M(1, 2); - const number2 t111 = M(3, 0) * M(1, 3); - const number2 t115 = M(0, 0) * M(2, 2); - const number2 t117 = M(0, 0) * M(2, 3); - const number2 t119 = M(2, 0) * M(0, 2); - const number2 t121 = M(2, 0) * M(0, 3); - const number2 t123 = M(3, 0) * M(0, 2); - const number2 t125 = M(3, 0) * M(0, 3); - const number2 t129 = M(0, 0) * M(1, 2); - const number2 t131 = M(0, 0) * M(1, 3); - const number2 t133 = M(1, 0) * M(0, 2); - const number2 t135 = M(1, 0) * M(0, 3); - (*this)(0, 0) = - (M(1, 1) * M(2, 2) * M(3, 3) - M(1, 1) * M(2, 3) * M(3, 2) - - M(2, 1) * M(1, 2) * M(3, 3) + M(2, 1) * M(1, 3) * M(3, 2) + - M(3, 1) * M(1, 2) * M(2, 3) - M(3, 1) * M(1, 3) * M(2, 2)) * - t65; - (*this)(0, 1) = - -(M(0, 1) * M(2, 2) * M(3, 3) - M(0, 1) * M(2, 3) * M(3, 2) - - t71 * M(3, 3) + t73 * M(3, 2) + t75 * M(2, 3) - t77 * M(2, 2)) * - t65; - (*this)(0, 2) = (t81 * M(3, 3) - t83 * M(3, 2) - t85 * M(3, 3) + - t87 * M(3, 2) + t75 * M(1, 3) - t77 * M(1, 2)) * - t65; - (*this)(0, 3) = -(t81 * M(2, 3) - t83 * M(2, 2) - t85 * M(2, 3) + - t87 * M(2, 2) + t71 * M(1, 3) - t73 * M(1, 2)) * - t65; - (*this)(1, 0) = - -(t101 * M(3, 3) - t103 * M(3, 2) - t105 * M(3, 3) + - t107 * M(3, 2) + t109 * M(2, 3) - t111 * M(2, 2)) * - t65; - (*this)(1, 1) = (t115 * M(3, 3) - t117 * M(3, 2) - t119 * M(3, 3) + - t121 * M(3, 2) + t123 * M(2, 3) - t125 * M(2, 2)) * - t65; - (*this)(1, 2) = - -(t129 * M(3, 3) - t131 * M(3, 2) - t133 * M(3, 3) + - t135 * M(3, 2) + t123 * M(1, 3) - t125 * M(1, 2)) * - t65; - (*this)(1, 3) = (t129 * M(2, 3) - t131 * M(2, 2) - t133 * M(2, 3) + - t135 * M(2, 2) + t119 * M(1, 3) - t121 * M(1, 2)) * - t65; - (*this)(2, 0) = (t32 * M(3, 3) - t103 * M(3, 1) - t46 * M(3, 3) + - t107 * M(3, 1) + t57 * M(2, 3) - t111 * M(2, 1)) * - t65; - (*this)(2, 1) = -(t19 * M(3, 3) - t117 * M(3, 1) - t43 * M(3, 3) + - t121 * M(3, 1) + t54 * M(2, 3) - t125 * M(2, 1)) * - t65; - (*this)(2, 2) = (t14 * M(3, 3) - t131 * M(3, 1) - t29 * M(3, 3) + - t135 * M(3, 1) + t54 * M(1, 3) - t125 * M(1, 1)) * - t65; - (*this)(2, 3) = -(t14 * M(2, 3) - t131 * M(2, 1) - t29 * M(2, 3) + - t135 * M(2, 1) + t43 * M(1, 3) - t121 * M(1, 1)) * - t65; - (*this)(3, 0) = -(t32 * M(3, 2) - t101 * M(3, 1) - t46 * M(3, 2) + - t105 * M(3, 1) + t57 * M(2, 2) - t109 * M(2, 1)) * - t65; - (*this)(3, 1) = (t19 * M(3, 2) - t115 * M(3, 1) - t43 * M(3, 2) + - t119 * M(3, 1) + t54 * M(2, 2) - t123 * M(2, 1)) * - t65; - (*this)(3, 2) = -(t14 * M(3, 2) - t129 * M(3, 1) - t29 * M(3, 2) + - t133 * M(3, 1) + t54 * M(1, 2) - t123 * M(1, 1)) * - t65; - (*this)(3, 3) = (t14 * M(2, 2) - t129 * M(2, 1) - t29 * M(2, 2) + - t133 * M(2, 1) + t43 * M(1, 2) - t119 * M(1, 1)) * - t65; - - break; - } - - - default: - // if no inversion is - // hardcoded, fall back - // to use the - // Gauss-Jordan algorithm + }; + + case 4: + { + // Initially derived from the following maple script + // + // with (linalg); + // a:=matrix(4,4); + // evalm(a); + // ai:=inverse(a); + // readlib(C); + // C(ai,optimized,filename=x4); + // + // but then combined re-occurring terms via distributive law in an + // FMA-friendly format + const number2 M00 = M(0, 0); + const number2 M01 = M(0, 1); + const number2 M02 = M(0, 2); + const number2 M03 = M(0, 3); + const number2 M10 = M(1, 0); + const number2 M11 = M(1, 1); + const number2 M12 = M(1, 2); + const number2 M13 = M(1, 3); + const number2 M20 = M(2, 0); + const number2 M21 = M(2, 1); + const number2 M22 = M(2, 2); + const number2 M23 = M(2, 3); + const number2 M30 = M(3, 0); + const number2 M31 = M(3, 1); + const number2 M32 = M(3, 2); + const number2 M33 = M(3, 3); + + const number2 t14 = M00 * M11 - M10 * M01; + const number2 t15 = M22 * M33 - M23 * M32; + const number2 t19 = M00 * M21 - M20 * M01; + const number2 t20 = M12 * M33 - M13 * M32; + const number2 t24 = M00 * M31 - M30 * M01; + const number2 t25 = M12 * M23 - M13 * M22; + const number2 t32 = M10 * M21 - M20 * M11; + const number2 t33 = M02 * M33 - M03 * M32; + const number2 t37 = M10 * M31 - M30 * M11; + const number2 t38 = M02 * M23 - M03 * M22; + const number2 t49 = M20 * M31 - M30 * M21; + const number2 t50 = M02 * M13 - M03 * M12; + const number2 det = t14 * t15 - t19 * t20 + t24 * t25 + t32 * t33 - + t37 * t38 + t49 * t50; + const number2 inv_det = number2(1.0) / det; + const number2 t81 = M01 * M12 - M02 * M11; + const number2 t83 = M01 * M13 - M03 * M11; + const number2 t93 = M01 * M22 - M02 * M21; + const number2 t95 = M11 * M23 - M13 * M21; + const number2 t97 = M01 * M23 - M03 * M21; + const number2 t99 = M11 * M22 - M12 * M21; + const number2 t101 = M10 * M22 - M20 * M12; + const number2 t103 = M10 * M23 - M20 * M13; + const number2 t115 = M00 * M22 - M20 * M02; + const number2 t117 = M00 * M23 - M20 * M03; + const number2 t129 = M00 * M12 - M10 * M02; + const number2 t131 = M00 * M13 - M10 * M03; + + (*this)(0, 0) = (M11 * t15 - M21 * t20 + M31 * t25) * inv_det; + (*this)(0, 1) = -(M01 * t15 - M21 * t33 + M31 * t38) * inv_det; + (*this)(0, 2) = (M33 * t81 - M32 * t83 + M31 * t50) * inv_det; + (*this)(0, 3) = -(M23 * t81 - M22 * t83 + M21 * t50) * inv_det; + (*this)(1, 0) = -(M33 * t101 - M32 * t103 + M30 * t25) * inv_det; + (*this)(1, 1) = (M33 * t115 - M32 * t117 + M30 * t38) * inv_det; + (*this)(1, 2) = -(M33 * t129 - M32 * t131 + M30 * t50) * inv_det; + (*this)(1, 3) = (M23 * t129 - M22 * t131 + M20 * t50) * inv_det; + (*this)(2, 0) = (M33 * t32 - M31 * t103 + M30 * t95) * inv_det; + (*this)(2, 1) = -(M33 * t19 - M31 * t117 + M30 * t97) * inv_det; + (*this)(2, 2) = (M33 * t14 - M31 * t131 + M30 * t83) * inv_det; + (*this)(2, 3) = -(M23 * t14 - M21 * t131 + M20 * t83) * inv_det; + (*this)(3, 0) = -(M32 * t32 - M31 * t101 + M30 * t99) * inv_det; + (*this)(3, 1) = (M32 * t19 - M31 * t115 + M30 * t93) * inv_det; + (*this)(3, 2) = -(M32 * t14 - M31 * t129 + M30 * t81) * inv_det; + (*this)(3, 3) = (M22 * t14 - M21 * t129 + M20 * t81) * inv_det; + + break; + } + + + default: + // if no inversion is + // hardcoded, fall back + // to use the + // Gauss-Jordan algorithm + if (!PointerComparison::equal(&M, this)) *this = M; - gauss_jordan(); - }; + gauss_jordan(); + }; }