[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sun, 6 May 2018 14:55:19 -0400 (EDT) |
branch: devel-yves-generic-assembly-modifs
commit 716a4e92320bcaea60e710ddb0e86d155798a57f
Author: Yves Renard <address@hidden>
Date: Sun May 6 12:18:23 2018 +0200
Adding operator Contract to the assembly language : some specific tensor
contractions
---
doc/sphinx/source/userdoc/gasm_high.rst | 14 +-
interface/tests/python/Makefile.am | 2 +
interface/tests/python/check_asm.py | 62 ++
src/getfem/getfem_generic_assembly_tree.h | 1 +
src/getfem_generic_assembly_compile_and_exec.cc | 762 +++++++++++++------
...fem_generic_assembly_functions_and_operators.cc | 3 +-
src/getfem_generic_assembly_interpolation.cc | 4 +-
src/getfem_generic_assembly_semantic.cc | 803 +++++++++++++--------
src/getfem_generic_assembly_tree.cc | 11 +-
9 files changed, 1099 insertions(+), 563 deletions(-)
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst
b/doc/sphinx/source/userdoc/gasm_high.rst
index 34448ca..c879f7b 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -58,9 +58,9 @@ A specific language has been developed to describe the weak
formulation of bound
- ``Reshape(t, i, j, ...)``: Reshape a vector/matrix/tensor. Note that all
tensors in |gf| are stored in the Fortran order.
- - A certain number of linear and nonlinear operators (``Trace``, ``Norm``,
``Det``, ``Deviator``, ...). The nonlinear operators cannot be applied to test
functions.
+ - A certain number of linear and nonlinear operators (``Trace``, ``Norm``,
``Det``, ``Deviator``, ``Contract``, ...). The nonlinear operators cannot be
applied to test functions.
- - ``Diff(expression, variable)``: The possibility to explicity differentiate
an expression with respect to a variable
+ - ``Diff(expression, variable)``: The possibility to explicity differentiate
an expression with respect to a variable.
- Possiblility of macro definition (in the model, the ga_workspace object or
directly in the assembly string). The macros should be some valid expressions
that are expanded inline at the lexical analysis phase (if they are used
several times, the computation is automatically factorized at the compilation
stage).
@@ -460,9 +460,9 @@ A certain number of binary operations between tensors are
available:
- ``/`` stands for the division by a scalar.
- - ``.`` stands for the scalar product of vectors, or more generally to the
reduction of a tensor with respect to the last index with a vector. Note that
``*`` and ``.`` are equivalent for matrix-vector multiplication.
+ - ``.`` stands for the scalar product of vectors, or more generally to the
contraction of a tensor with respect to its last index with a vector or with
the first index of another tensor. Note that ``*`` and ``.`` are equivalent for
matrix-vector or matrix-matrix multiplication.
- - ``:`` stands for the the |Frobenius| product of matrices or more
generally to the reduction of a tensor with respect to the two last indices
with a matrix. Note that ``*`` and ``:`` are equivalent for (fourth order
tensor)-matrix multiplication.
+ - ``:`` stands for the the |Frobenius| product of matrices or more
generally to the contraction of a tensor with respect to the two last indices
with a matrix. Note that ``*`` and ``:`` are equivalent for (fourth order
tensor)-matrix multiplication.
- ``.*`` stands for the multiplication of two vectors/matrix/tensor
componentwise.
@@ -470,6 +470,10 @@ A certain number of binary operations between tensors are
available:
- address@hidden stands for the tensor product.
+ - ``Contract(A, i, B, j)`` stands for the contraction of tensors A and B
with respect to the ith index of A and jth index of B. The first index is
numbered 1. For instance ``Contract(V,1,W,1)`` is equivalent to ``V.W`` for two
vectors ``V`` and ``W``.
+
+ - ``Contract(A, i, j, B, k, l)`` stands for the double contraction of
tensors A and B with respect to indices i,j of A and indices k,l of B. The
first index is numbered 1. For instance ``Contract(A,1,2,B,1,2)`` is equivalent
to ``A:B`` for two matrices ``A`` and ``B``.
+
Unary operators
---------------
@@ -478,7 +482,9 @@ Unary operators
- ``'`` stands for the transpose of a matrix or line view of a vector.
+ - ``Contract(A, i, j)`` stands for the contraction of tensor A with respect
to its ith and jth indices. The first index is numbered 1. For instance,
``Contract(A, 1, 2)`` is equivalent to ``Trace(A)`` for a matrix ``A``.
+
Parentheses
-----------
diff --git a/interface/tests/python/Makefile.am
b/interface/tests/python/Makefile.am
index 40bac60..64c4b46 100644
--- a/interface/tests/python/Makefile.am
+++ b/interface/tests/python/Makefile.am
@@ -26,6 +26,7 @@ EXTRA_DIST= \
check_export.py \
check_global_functions.py \
check_levelset.py \
+ check_asm.py \
demo_crack.py \
demo_fictitious_domains.py \
demo_laplacian.py \
@@ -54,6 +55,7 @@ if BUILDPYTHON
TESTS = \
check_export.py \
check_global_functions.py \
+ check_asm.py \
demo_wave.py \
demo_laplacian.py \
$(optpy)
diff --git a/interface/tests/python/check_asm.py
b/interface/tests/python/check_asm.py
new file mode 100644
index 0000000..ad8bf86
--- /dev/null
+++ b/interface/tests/python/check_asm.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2018-2018 Yves Renard.
+#
+# This file is a part of GetFEM++
+#
+# GetFEM++ is free software; you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation; either version 2.1 of the License, or
+# (at your option) any later version.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+# License for more details.
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+############################################################################
+""" test high generic assembly language.
+
+ This program is used to check that Python-GetFEM interface, and more
+ generally GetFEM are working. It focuses on testing some operations
+ of the high generic assembly language.
+
+ $Id$
+"""
+import numpy as np
+import getfem as gf
+import os
+
+
+NX = 4
+m = gf.Mesh('triangles grid', np.arange(0,1+1./NX,1./NX),
+ np.arange(0,1+1./NX,1./NX)) # Structured mesh
+fem = gf.Fem('FEM_PK(2,1)')
+mfu = gf.MeshFem(m, 1); mfu.set_fem(fem) # Lagrange P1 scalar fem
+mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
+mfv = gf.MeshFem(m, 3); mfv.set_fem(fem) # Lagrange P1 vector fem
+
+
+U = mfu.eval('x+y')
+V = mfv.eval('[x*y, x*y, x*y]')
+
+
+md = gf.Model('real')
+
+md.add_fem_variable('u', mfu)
+md.set_variable('u', U)
+md.add_fem_variable('v', mfv)
+md.set_variable('v', V)
+
+
+print gf.asm('generic', mim, 0, "Def P(a):=a*(a'); Print(P(Grad_v))", -1, md)
+print gf.asm('generic', mim, 1, "Print(Grad_Test_u)", -1, md)
+print gf.asm('generic', mim, 0, "Def P(a):=a*(a'); Contract(P(Grad_v), 1, 2)",
-1, md)
+
+
+
+
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index 2d05fdc..34a5cde 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -131,6 +131,7 @@ namespace getfem {
GA_NODE_MACRO_PARAM,
GA_NODE_PARAMS,
GA_NODE_RESHAPE,
+ GA_NODE_CONTRACT,
GA_NODE_ALLINDICES,
GA_NODE_C_MATRIX,
GA_NODE_X,
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index d7b14dc..9f2c794 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -1929,8 +1929,7 @@ namespace getfem {
struct ga_instruction_dotmult_spec : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: specific componentwise "
- "multiplication");
+ GA_DEBUG_INFO("Instruction: specific componentwise multiplication");
size_type s2_1 = tc2.sizes()[0], s2_2 = tc2.size() / s2_1;
size_type s1_1 = tc1.size() / s2_2;
@@ -1947,69 +1946,290 @@ namespace getfem {
: t(t_), tc1(tc1_), tc2(tc2_) {}
};
- // Performs Amij Bjk -> Cmik. To be optimized
+ // Performs Amijik -> Cmjk. To be optimized
+ struct ga_instruction_contract_1_1 : public ga_instruction {
+ base_tensor &t, &tc1;
+ size_type nn, ii2, ii3;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: single contraction on a single tensor");
+
+ size_type ii1 = tc1.size() / (nn*ii2*ii3);
+
+ base_tensor::iterator it = t.begin();
+ for (size_type i = 0; i < ii3; ++i)
+ for (size_type j = 0; j < ii2; ++j)
+ for (size_type k = 0; k < ii1; ++k, ++it) {
+ *it = scalar_type(0);
+ size_type pre_ind = k+j*ii1*nn+i*ii1*nn*ii2*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += tc1[pre_ind+n*ii1+n*ii1*nn*ii2];
+ }
+
+ GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
+ return 0;
+ }
+ ga_instruction_contract_1_1(base_tensor &t_, base_tensor &tc1_,
+ size_type n_, size_type i2_, size_type i3_)
+ : t(t_), tc1(tc1_), nn(n_), ii2(i2_), ii3(i3_) {}
+ };
+
+ // Performs Amijk Bnljp -> Cmniklp. To be optimized
+ struct ga_instruction_contract_2_1 : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ size_type nn, ii1, ii2, ii3, ii4;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: single contraction on two tensors");
+
+ size_type ift1 = tc1.size() / (nn*ii1*ii2);
+ size_type ift2 = tc2.size() / (nn*ii3*ii4);
+
+ base_tensor::iterator it = t.begin();
+ for (size_type i = 0; i < ii4; ++i)
+ for (size_type j = 0; j < ii3; ++j)
+ for (size_type k = 0; k < ii2; ++k)
+ for (size_type l = 0; l < ii1; ++l)
+ for (size_type p = 0; p < ift2; ++p)
+ for (size_type q = 0; q < ift1; ++q, ++it) {
+ *it = scalar_type(0);
+ size_type ind1 = q+l*ift1+k*ift1*ii1*nn;
+ size_type ind2 = p+j*ift2+i*ift2*ii3*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += tc1[ind1+n*ift1*ii1] * tc2[ind2+n*ift2*ii3];
+ }
+
+ GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
+ return 0;
+ }
+ ga_instruction_contract_2_1(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_,
+ size_type n_, size_type i1_, size_type i2_,
+ size_type i3_, size_type i4_)
+ : t(t_), tc1(tc1_), tc2(tc2_), nn(n_),
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_) {}
+ };
+
+ // Performs Amijk Bnljp -> Cnmiklp. To be optimized
+ struct ga_instruction_contract_2_1_rev : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ size_type nn, ii1, ii2, ii3, ii4;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: single contraction on two tensors");
+
+ size_type ift1 = tc1.size() / (nn*ii1*ii2);
+ size_type ift2 = tc2.size() / (nn*ii3*ii4);
+
+ base_tensor::iterator it = t.begin();
+ for (size_type i = 0; i < ii4; ++i)
+ for (size_type j = 0; j < ii3; ++j)
+ for (size_type k = 0; k < ii2; ++k)
+ for (size_type l = 0; l < ii1; ++l)
+ for (size_type q = 0; q < ift1; ++q)
+ for (size_type p = 0; p < ift2; ++p, ++it) {
+ *it = scalar_type(0);
+ size_type ind1 = q+l*ift1+k*ift1*ii1*nn;
+ size_type ind2 = p+j*ift2+i*ift2*ii3*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += tc1[ind1+n*ift1*ii1] * tc2[ind2+n*ift2*ii3];
+ }
+
+ GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
+ return 0;
+ }
+ ga_instruction_contract_2_1_rev(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_,
+ size_type n_, size_type i1_, size_type i2_,
+ size_type i3_, size_type i4_)
+ : t(t_), tc1(tc1_), tc2(tc2_), nn(n_),
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_) {}
+ };
+
+ // Performs Amijklp Bnqjrls -> Cmnikpqrs. To be optimized
+ struct ga_instruction_contract_2_2 : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ size_type nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6;
+ bool inv_tc2;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: single contraction on two tensors");
+
+ size_type ift1 = tc1.size() / (nn1*nn2*ii1*ii2*ii3);
+ size_type ift2 = tc2.size() / (nn1*nn2*ii3*ii4*ii5);
+
+ size_type sn1 = ift2*ii4, sn2 = ift2*ii4*nn1*ii5;
+ if (inv_tc2) std::swap(sn1, sn2);
+
+ base_tensor::iterator it = t.begin();
+ for (size_type i = 0; i < ii6; ++i)
+ for (size_type j = 0; j < ii5; ++j)
+ for (size_type k = 0; k < ii4; ++k)
+ for (size_type l = 0; l < ii3; ++l)
+ for (size_type p = 0; p < ii2; ++p)
+ for (size_type q = 0; q < ii1; ++q)
+ for (size_type r = 0; r < ift2; ++r)
+ for (size_type s = 0; s < ift1; ++s, ++it) {
+ *it = scalar_type(0);
+ size_type ind1
+ = s+q*ift1+p*ift1*ii1*nn1+l*ift1*ii1*nn1*ii2*nn2;
+ size_type ind2
+ = r+k*ift2+j*ift2*ii4*nn1+i*ift2*ii4*nn1*ii5*nn2;
+ for (size_type n1 = 0; n1 < nn1; ++n1)
+ for (size_type n2 = 0; n2 < nn2; ++n2)
+ *it += tc1[ind1+n1*ift1*ii1+n2*ift1*ii1*nn1*ii2]
+ * tc2[ind2+n1*sn1+n2*sn2];
+ }
+
+ GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
+ return 0;
+ }
+ ga_instruction_contract_2_2(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_,
+ size_type n1_, size_type n2_,
+ size_type i1_, size_type i2_, size_type i3_,
+ size_type i4_, size_type i5_, size_type i6_,
+ bool intc2)
+ : t(t_), tc1(tc1_), tc2(tc2_), nn1(n1_), nn2(n2_),
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_), ii5(i5_), ii6(i6_),
+ inv_tc2(intc2) {}
+ };
+
+ // Performs Amijklp Bnqjrls -> Cnmikpqrs. To be optimized
+ struct ga_instruction_contract_2_2_rev : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ size_type nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6;
+ bool inv_tc2;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: single contraction on two tensors");
+
+ size_type ift1 = tc1.size() / (nn1*nn2*ii1*ii2*ii3);
+ size_type ift2 = tc2.size() / (nn1*nn2*ii3*ii4*ii5);
+
+ size_type sn1 = ift2*ii4, sn2 = ift2*ii4*nn1*ii5;
+ if (inv_tc2) std::swap(sn1, sn2);
+
+ base_tensor::iterator it = t.begin();
+ for (size_type i = 0; i < ii6; ++i)
+ for (size_type j = 0; j < ii5; ++j)
+ for (size_type k = 0; k < ii4; ++k)
+ for (size_type l = 0; l < ii3; ++l)
+ for (size_type p = 0; p < ii2; ++p)
+ for (size_type q = 0; q < ii1; ++q)
+ for (size_type s = 0; s < ift1; ++s)
+ for (size_type r = 0; r < ift2; ++r, ++it) {
+ *it = scalar_type(0);
+ size_type ind1
+ = s+q*ift1+p*ift1*ii1*nn1+l*ift1*ii1*nn1*ii2*nn2;
+ size_type ind2
+ = r+k*ift2+j*ift2*ii4*nn1+i*ift2*ii4*nn1*ii5*nn2;
+ for (size_type n1 = 0; n1 < nn1; ++n1)
+ for (size_type n2 = 0; n2 < nn2; ++n2)
+ *it += tc1[ind1+n1*ift1*ii1+n2*ift1*ii1*nn1*ii2]
+ * tc2[ind2+n1*sn1+n2*sn2];
+ }
+
+ GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
+ return 0;
+ }
+ ga_instruction_contract_2_2_rev(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_,
+ size_type n1_, size_type n2_,
+ size_type i1_, size_type i2_, size_type i3_,
+ size_type i4_, size_type i5_, size_type i6_,
+ bool intc2)
+ : t(t_), tc1(tc1_), tc2(tc2_), nn1(n1_), nn2(n2_),
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_), ii5(i5_), ii6(i6_),
+ inv_tc2(intc2) {}
+ };
+
+
+ // Performs Amj Bjk -> Cmk. To be optimized
struct ga_instruction_matrix_mult : public ga_instruction {
base_tensor &t, &tc1, &tc2;
+ size_type n;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: matrix multiplication");
- size_type order = tc2.sizes().size();
- size_type s2_1 = tc2.sizes()[order-2];
- size_type s2_2 = tc2.sizes()[order-1];
- size_type s1 = tc1.size() / s2_1;
- size_type s2 = tc2.size() / (s2_1*s2_2);
+ GA_DEBUG_INFO("Instruction: order one contraction "
+ "(dot product or matrix multiplication)");
+
+ size_type s1 = tc1.size() / n;
+ size_type s2 = tc2.size() / n;
base_tensor::iterator it = t.begin();
- for (size_type k = 0; k < s2_2; ++k)
- for (size_type i = 0; i < s1; ++i)
- for (size_type m = 0; m < s2; ++m, ++it) {
+ for (size_type k = 0; k < s2; ++k)
+ for (size_type i = 0; i < s1; ++i, ++it) {
*it = scalar_type(0);
- for (size_type j = 0; j < s2_1; ++j)
- *it += tc1[i+j*s1] * tc2[m+j*s2+k*s2_1*s2];
+ for (size_type j = 0; j < n; ++j)
+ *it += tc1[i+j*s1] * tc2[j+k*n];
}
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_matrix_mult(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_)
- : t(t_), tc1(tc1_), tc2(tc2_) {}
+ base_tensor &tc2_, size_type n_)
+ : t(t_), tc1(tc1_), tc2(tc2_), n(n_) {}
};
// Performs Amij Bnjk -> Cmnik. To be optimized
struct ga_instruction_matrix_mult_spec : public ga_instruction {
base_tensor &t, &tc1, &tc2;
+ size_type n, m, p; // tc1 of size q*m*n, tc2 of size l*n*p
+ // t of size q*l*m*p
virtual int exec() {
- GA_DEBUG_INFO("Instruction: specific matrix multiplication");
- size_type s1_1 = tc1.sizes()[0];
- size_type s1_2 = tc1.sizes()[1];
- // size_type s1_3 = tc1.sizes()[2];
- size_type s2_1 = tc2.sizes()[0];
- size_type s2_2 = tc2.sizes()[1];
- size_type s2_3 = tc2.sizes()[2];
-
+ GA_DEBUG_INFO("Instruction: specific order one contraction "
+ "(dot product or matrix multiplication)");
+ size_type q = tc1.size() / (m * n);
+ size_type l = tc2.size() / (p * n);
+
base_tensor::iterator it = t.begin();
- for (size_type k = 0; k < s2_3; ++k)
- for (size_type i = 0; i < s1_2; ++i)
- for (size_type n = 0; n < s2_1; ++n)
- for (size_type m = 0; m < s1_1; ++m, ++it) {
- *it = scalar_type(0);
- for (size_type j = 0; j < s2_2; ++j)
- *it += tc1[m+i*s1_1+j*s1_1*s1_2] * tc2[n+j*s2_1+k*s2_1*s2_2];
- }
+ for (size_type r = 0; r < p; ++r)
+ for (size_type k = 0; k < m; ++k)
+ for (size_type j = 0; j < l; ++j)
+ for (size_type i = 0; i < q; ++i, ++it) {
+ *it = scalar_type(0);
+ for (size_type s = 0; s < n; ++s)
+ *it += tc1[i+k*q+s*q*m] * tc2[j+s*l+r*l*n];
+ }
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_matrix_mult_spec(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_)
- : t(t_), tc1(tc1_), tc2(tc2_) {}
+ base_tensor &tc2_, size_type n_,
+ size_type m_, size_type p_)
+ : t(t_), tc1(tc1_), tc2(tc2_), n(n_), m(m_), p(p_) {}
};
+ // Performs Amij Bnjk -> Cnmik. To be optimized
+ struct ga_instruction_matrix_mult_spec2 : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ size_type n, m, p; // tc1 of size q*m*n, tc2 of size l*n*p
+ // t of size l*q*m*p
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: specific order one contraction "
+ "(dot product or matrix multiplication)");
+ size_type q = tc1.size() / (m * n);
+ size_type l = tc2.size() / (p * n);
+
+ base_tensor::iterator it = t.begin();
+ for (size_type r = 0; r < p; ++r)
+ for (size_type k = 0; k < m; ++k)
+ for (size_type i = 0; i < q; ++i)
+ for (size_type j = 0; j < l; ++j, ++it) {
+ *it = scalar_type(0);
+ for (size_type s = 0; s < n; ++s)
+ *it += tc1[i+k*q+s*q*m] * tc2[j+s*l+r*l*n];
+ }
+ GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
+ return 0;
+ }
+ ga_instruction_matrix_mult_spec2(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_, size_type n_,
+ size_type m_, size_type p_)
+ : t(t_), tc1(tc1_), tc2(tc2_), n(n_), m(m_), p(p_) {}
+ };
// Performs Ani Bmi -> Cmn
- struct ga_instruction_reduction : public ga_instruction {
+ struct ga_instruction_contraction : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type nn;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: reduction operation of size " << nn);
+ GA_DEBUG_INFO("Instruction: contraction operation of size " << nn);
#if GA_USES_BLAS
long m = int(tc1.size()/nn), k = int(nn), n = int(tc2.size()/nn);
long lda = m, ldb = n, ldc = m;
@@ -2040,17 +2260,17 @@ namespace getfem {
#endif
return 0;
}
- ga_instruction_reduction(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_contraction(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_) {}
};
// Performs Ani Bmi -> Cmn
- struct ga_instruction_reduction_opt0_2 : public ga_instruction {
+ struct ga_instruction_contraction_opt0_2 : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type n, q;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: reduction operation of size " << n*q <<
+ GA_DEBUG_INFO("Instruction: contraction operation of size " << n*q <<
" optimized for vectorized second tensor of type 2");
size_type nn = n*q, s1 = tc1.size()/nn, s2 = tc2.size()/nn, s2_q = s2/q;
size_type s1_qq = s1*q, s2_qq = s2*q;
@@ -2073,12 +2293,12 @@ namespace getfem {
}
}
// base_tensor u = t;
- // ga_instruction_reduction toto(t, tc1, tc2, n*q);
+ // ga_instruction_contraction toto(t, tc1, tc2, n*q);
// toto.exec();
// GMM_ASSERT1(gmm::vect_dist2(t.as_vector(), u.as_vector()) < 1E-9,
"Erroneous");
return 0;
}
- ga_instruction_reduction_opt0_2(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_contraction_opt0_2(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_,
size_type q_)
: t(t_), tc1(tc1_), tc2(tc2_), n(n_), q(q_) {}
@@ -2086,11 +2306,11 @@ namespace getfem {
// Performs Ani Bmi -> Cmn
template <int N>
- struct ga_instruction_reduction_opt0_2_unrolled : public ga_instruction {
+ struct ga_instruction_contraction_opt0_2_unrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type q;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N*q
+ GA_DEBUG_INFO("Instruction: unrolled contraction operation of size " <<
N*q
<< " optimized for vectorized second tensor of type 2");
size_type nn = N*q, s1 = tc1.size()/nn, s2 = tc2.size()/nn, s2_q = s2/q;
size_type s1_qq = s1*q, s2_qq = s2*q;
@@ -2114,17 +2334,17 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt0_2_unrolled(base_tensor &t_, base_tensor
&tc1_,
+ ga_instruction_contraction_opt0_2_unrolled(base_tensor &t_, base_tensor
&tc1_,
base_tensor &tc2_, size_type q_)
: t(t_), tc1(tc1_), tc2(tc2_), q(q_) {}
};
// Performs Ani Bmi -> Cmn
template <int N, int Q>
- struct ga_instruction_reduction_opt0_2_dunrolled : public ga_instruction {
+ struct ga_instruction_contraction_opt0_2_dunrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N*Q
+ GA_DEBUG_INFO("Instruction: unrolled contraction operation of size " <<
N*Q
<< " optimized for vectorized second tensor of type 2");
size_type s1 = tc1.size()/(N*Q), s2 = tc2.size()/(N*Q), s2_q = s2/Q;
size_type s1_qq = s1*Q, s2_qq = s2*Q;
@@ -2148,17 +2368,17 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt0_2_dunrolled
+ ga_instruction_contraction_opt0_2_dunrolled
(base_tensor &t_, base_tensor &tc1_, base_tensor &tc2_)
: t(t_), tc1(tc1_), tc2(tc2_) {}
};
// Performs Ani Bmi -> Cmn
- struct ga_instruction_reduction_opt2_0 : public ga_instruction {
+ struct ga_instruction_contraction_opt2_0 : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type n, q;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: reduction operation of size " << n*q <<
+ GA_DEBUG_INFO("Instruction: contraction operation of size " << n*q <<
" optimized for vectorized second tensor of type 2");
size_type nn = n*q, s1 = tc1.size()/nn, s2 = tc2.size()/nn;
size_type s1_q = s1/q, s1_qq = s1*q, s2_qq = s2*q;
@@ -2180,7 +2400,7 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt2_0(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_contraction_opt2_0(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_,
size_type q_)
: t(t_), tc1(tc1_), tc2(tc2_), n(n_), q(q_) { }
@@ -2188,11 +2408,11 @@ namespace getfem {
// Performs Ani Bmi -> Cmn
template <int N>
- struct ga_instruction_reduction_opt2_0_unrolled : public ga_instruction {
+ struct ga_instruction_contraction_opt2_0_unrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type q;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N*q
+ GA_DEBUG_INFO("Instruction: unrolled contraction operation of size " <<
N*q
<< " optimized for vectorized second tensor of type 2");
size_type nn = N*q, s1 = tc1.size()/nn, s2 = tc2.size()/nn;
size_type s1_q = s1/q, s1_qq = s1*q, s2_qq = s2*q;
@@ -2213,17 +2433,17 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt2_0_unrolled(base_tensor &t_, base_tensor
&tc1_,
+ ga_instruction_contraction_opt2_0_unrolled(base_tensor &t_, base_tensor
&tc1_,
base_tensor &tc2_, size_type q_)
: t(t_), tc1(tc1_), tc2(tc2_), q(q_) {}
};
// Performs Ani Bmi -> Cmn
template <int N, int Q>
- struct ga_instruction_reduction_opt2_0_dunrolled : public ga_instruction {
+ struct ga_instruction_contraction_opt2_0_dunrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N*Q
+ GA_DEBUG_INFO("Instruction: unrolled contraction operation of size " <<
N*Q
<< " optimized for vectorized second tensor of type 2");
size_type s1 = tc1.size()/(N*Q), s2 = tc2.size()/(N*Q);
size_type s1_q = s1/Q, s1_qq = s1*Q, s2_qq = s2*Q;
@@ -2244,17 +2464,17 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt2_0_dunrolled
+ ga_instruction_contraction_opt2_0_dunrolled
(base_tensor &t_, base_tensor &tc1_, base_tensor &tc2_)
: t(t_), tc1(tc1_), tc2(tc2_) {}
};
// Performs Ani Bmi -> Cmn
- struct ga_instruction_reduction_opt0_1 : public ga_instruction {
+ struct ga_instruction_contraction_opt0_1 : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type nn;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: reduction operation of size " << nn <<
+ GA_DEBUG_INFO("Instruction: contraction operation of size " << nn <<
" optimized for vectorized second tensor of type 1");
size_type ss1=tc1.size(), s1 = ss1/nn, s2=tc2.size()/nn, s2_n=s2/nn;
@@ -2271,7 +2491,7 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt0_1(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_contraction_opt0_1(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_) {}
};
@@ -2289,10 +2509,10 @@ namespace getfem {
// Performs Ani Bmi -> Cmn
template <int N>
- struct ga_instruction_reduction_opt0_1_unrolled : public ga_instruction {
+ struct ga_instruction_contraction_opt0_1_unrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N
+ GA_DEBUG_INFO("Instruction: unrolled contraction operation of size " << N
<< " optimized for vectorized second tensor of type 1");
size_type s1 = tc1.size()/N, s2 = tc2.size()/N;
auto it = t.begin(), it1 = tc1.begin();
@@ -2303,17 +2523,17 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt0_1_unrolled(base_tensor &t_, base_tensor
&tc1_,
+ ga_instruction_contraction_opt0_1_unrolled(base_tensor &t_, base_tensor
&tc1_,
base_tensor &tc2_)
: t(t_), tc1(tc1_), tc2(tc2_) {}
};
// Performs Ani Bmi -> Cmn
- struct ga_instruction_reduction_opt1_1 : public ga_instruction {
+ struct ga_instruction_contraction_opt1_1 : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type nn;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: reduction operation of size " << nn <<
+ GA_DEBUG_INFO("Instruction: contraction operation of size " << nn <<
" optimized for both vectorized tensor of type 1");
size_type s1 = tc1.size()/nn, s2 = tc2.size()/nn, s2_1 = s2+1;
GA_DEBUG_ASSERT(t.size() == s2*s1, "Internal error");
@@ -2334,7 +2554,7 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_opt1_1(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_contraction_opt1_1(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_) {}
};
@@ -2353,11 +2573,11 @@ namespace getfem {
{ return (*it1)*(*it2); }
// Performs Ani Bmi -> Cmn. Unrolled operation.
- template<int N> struct ga_instruction_reduction_unrolled
+ template<int N> struct ga_instruction_contraction_unrolled
: public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N);
+ GA_DEBUG_INFO("Instruction: unrolled contraction operation of size " <<
N);
size_type s1 = tc1.size()/N, s2 = tc2.size()/N;
GA_DEBUG_ASSERT(t.size() == s1*s2, "Internal error, " << t.size()
<< " != " << s1 << "*" << s2);
@@ -2368,7 +2588,7 @@ namespace getfem {
}
return 0;
}
- ga_instruction_reduction_unrolled(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_contraction_unrolled(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_)
: t(t_), tc1(tc1_), tc2(tc2_) {}
};
@@ -2437,7 +2657,7 @@ namespace getfem {
: public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: doubly unrolled reduction operation of size "
+ GA_DEBUG_INFO("Instruction: doubly unrolled contraction operation of
size "
<< S2 << "x" << N);
size_type s1 = tc1.size()/N, s2 = tc2.size()/N;
GA_DEBUG_ASSERT(s2 == S2, "Internal error");
@@ -2456,7 +2676,7 @@ namespace getfem {
};
- pga_instruction ga_instruction_reduction_switch
+ pga_instruction ga_instruction_contraction_switch
(assembly_tensor &t_, assembly_tensor &tc1_, assembly_tensor &tc2_,
size_type n, bool &to_clear) {
base_tensor &t = t_.tensor(), &tc1 = tc1_.tensor(), &tc2 = tc2_.tensor();
@@ -2465,25 +2685,25 @@ namespace getfem {
tc1_.qdim() == n && tc2_.qdim() == n) {
to_clear = true;
t_.set_sparsity(10, tc1_.qdim());
- return std::make_shared<ga_instruction_reduction_opt1_1>(t, tc1, tc2, n);
+ return std::make_shared<ga_instruction_contraction_opt1_1>(t, tc1, tc2,
n);
}
if (tc2_.sparsity() == 1) {
switch(n) {
case 2:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<2>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<2>>
(t, tc1, tc2);
case 3:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<3>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<3>>
(t, tc1, tc2);
case 4:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<4>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<4>>
(t, tc1, tc2);
case 5:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<5>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<5>>
(t, tc1, tc2);
default:
- return std::make_shared<ga_instruction_reduction_opt0_1>(t,tc1,tc2, n);
+ return std::make_shared<ga_instruction_contraction_opt0_1>(t,tc1,tc2,
n);
}
}
if (tc2_.sparsity() == 2) {
@@ -2495,64 +2715,64 @@ namespace getfem {
switch (q2) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<1,2>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<1,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<1,3>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<1,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<1,4>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<1,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt0_2_unrolled<1>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<1>>
(t, tc1, tc2, q2);
}
case 2:
switch (q2) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<2,2>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<2,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<2,3>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<2,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<2,4>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<2,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt0_2_unrolled<2>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<2>>
(t, tc1, tc2, q2);
}
case 3:
switch (q2) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<3,2>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<3,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<3,3>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<3,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<3,4>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<3,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt0_2_unrolled<3>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<3>>
(t, tc1, tc2, q2);
}
case 4:
- return std::make_shared<ga_instruction_reduction_opt0_2_unrolled<4>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<4>>
(t, tc1, tc2, q2);
case 5:
- return std::make_shared<ga_instruction_reduction_opt0_2_unrolled<5>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<5>>
(t, tc1, tc2, q2);
default:
- return std::make_shared<ga_instruction_reduction_opt0_2>
+ return std::make_shared<ga_instruction_contraction_opt0_2>
(t,tc1,tc2,n2,q2);
}
}
@@ -2566,108 +2786,108 @@ namespace getfem {
switch (q1) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<1,2>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<1,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<1,3>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<1,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<1,4>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<1,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt2_0_unrolled<1>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<1>>
(t, tc1, tc2, q1);
}
case 2:
switch (q1) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<2,2>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<2,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<2,3>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<2,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<2,4>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<2,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt2_0_unrolled<2>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<2>>
(t, tc1, tc2, q1);
}
case 3:
switch (q1) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<3,2>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<3,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<3,3>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<3,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<3,4>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<3,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt2_0_unrolled<3>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<3>>
(t, tc1, tc2, q1);
}
- return std::make_shared<ga_instruction_reduction_opt2_0_unrolled<3>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<3>>
(t, tc1, tc2, q1);
case 4:
- return std::make_shared<ga_instruction_reduction_opt2_0_unrolled<4>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<4>>
(t, tc1, tc2, q1);
case 5:
- return std::make_shared<ga_instruction_reduction_opt2_0_unrolled<5>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<5>>
(t, tc1, tc2, q1);
default:
- return std::make_shared<ga_instruction_reduction_opt2_0>
+ return std::make_shared<ga_instruction_contraction_opt2_0>
(t,tc1,tc2, n1, q1);
}
}
}
switch(n) {
- case 2 : return std::make_shared<ga_instruction_reduction_unrolled< 2>>
+ case 2 : return std::make_shared<ga_instruction_contraction_unrolled< 2>>
(t, tc1, tc2);
- case 3 : return std::make_shared<ga_instruction_reduction_unrolled< 3>>
+ case 3 : return std::make_shared<ga_instruction_contraction_unrolled< 3>>
(t, tc1, tc2);
- case 4 : return std::make_shared<ga_instruction_reduction_unrolled< 4>>
+ case 4 : return std::make_shared<ga_instruction_contraction_unrolled< 4>>
(t, tc1, tc2);
- case 5 : return std::make_shared<ga_instruction_reduction_unrolled< 5>>
+ case 5 : return std::make_shared<ga_instruction_contraction_unrolled< 5>>
(t, tc1, tc2);
- case 6 : return std::make_shared<ga_instruction_reduction_unrolled< 6>>
+ case 6 : return std::make_shared<ga_instruction_contraction_unrolled< 6>>
(t, tc1, tc2);
- case 7 : return std::make_shared<ga_instruction_reduction_unrolled< 7>>
+ case 7 : return std::make_shared<ga_instruction_contraction_unrolled< 7>>
(t, tc1, tc2);
- case 8 : return std::make_shared<ga_instruction_reduction_unrolled< 8>>
+ case 8 : return std::make_shared<ga_instruction_contraction_unrolled< 8>>
(t, tc1, tc2);
- case 9 : return std::make_shared<ga_instruction_reduction_unrolled< 9>>
+ case 9 : return std::make_shared<ga_instruction_contraction_unrolled< 9>>
(t, tc1, tc2);
- case 10 : return std::make_shared<ga_instruction_reduction_unrolled<10>>
+ case 10 : return std::make_shared<ga_instruction_contraction_unrolled<10>>
(t, tc1, tc2);
- case 11 : return std::make_shared<ga_instruction_reduction_unrolled<11>>
+ case 11 : return std::make_shared<ga_instruction_contraction_unrolled<11>>
(t, tc1, tc2);
- case 12 : return std::make_shared<ga_instruction_reduction_unrolled<12>>
+ case 12 : return std::make_shared<ga_instruction_contraction_unrolled<12>>
(t, tc1, tc2);
- case 13 : return std::make_shared<ga_instruction_reduction_unrolled<13>>
+ case 13 : return std::make_shared<ga_instruction_contraction_unrolled<13>>
(t, tc1, tc2);
- case 14 : return std::make_shared<ga_instruction_reduction_unrolled<14>>
+ case 14 : return std::make_shared<ga_instruction_contraction_unrolled<14>>
(t, tc1, tc2);
- case 15 : return std::make_shared<ga_instruction_reduction_unrolled<15>>
+ case 15 : return std::make_shared<ga_instruction_contraction_unrolled<15>>
(t, tc1, tc2);
- case 16 : return std::make_shared<ga_instruction_reduction_unrolled<16>>
+ case 16 : return std::make_shared<ga_instruction_contraction_unrolled<16>>
(t, tc1, tc2);
- default : return std::make_shared<ga_instruction_reduction>
+ default : return std::make_shared<ga_instruction_contraction>
(t, tc1, tc2, n);
}
}
- pga_instruction ga_uniform_instruction_reduction_switch
+ pga_instruction ga_uniform_instruction_contraction_switch
(assembly_tensor &t_, assembly_tensor &tc1_, assembly_tensor &tc2_,
size_type n, bool &to_clear) {
base_tensor &t = t_.tensor(), &tc1 = tc1_.tensor(), &tc2 = tc2_.tensor();
@@ -2676,24 +2896,24 @@ namespace getfem {
tc1_.qdim() == n && tc2_.qdim() == n) {
to_clear = true;
t_.set_sparsity(10, tc1_.qdim());
- return std::make_shared<ga_instruction_reduction_opt1_1>(t,tc1,tc2,n);
+ return std::make_shared<ga_instruction_contraction_opt1_1>(t,tc1,tc2,n);
}
if (tc2_.sparsity() == 1) {
switch(n) {
case 2:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<2>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<2>>
(t, tc1, tc2);
case 3:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<3>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<3>>
(t, tc1, tc2);
case 4:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<4>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<4>>
(t, tc1, tc2);
case 5:
- return std::make_shared<ga_instruction_reduction_opt0_1_unrolled<5>>
+ return std::make_shared<ga_instruction_contraction_opt0_1_unrolled<5>>
(t, tc1, tc2);
default:
- return std::make_shared<ga_instruction_reduction_opt0_1>(t,tc1,tc2, n);
+ return std::make_shared<ga_instruction_contraction_opt0_1>(t,tc1,tc2,
n);
}
}
if (tc2_.sparsity() == 2) {
@@ -2705,64 +2925,64 @@ namespace getfem {
switch (q2) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<1,2>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<1,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<1,3>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<1,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<1,4>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<1,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt0_2_unrolled<1>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<1>>
(t, tc1, tc2, q2);
}
case 2:
switch (q2) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<2,2>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<2,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<2,3>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<2,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<2,4>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<2,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt0_2_unrolled<2>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<2>>
(t, tc1, tc2, q2);
}
case 3:
switch (q2) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<3,2>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<3,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<3,3>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<3,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt0_2_dunrolled<3,4>>
+
std::make_shared<ga_instruction_contraction_opt0_2_dunrolled<3,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt0_2_unrolled<3>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<3>>
(t, tc1, tc2, q2);
}
case 4:
- return std::make_shared<ga_instruction_reduction_opt0_2_unrolled<4>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<4>>
(t, tc1, tc2, q2);
case 5:
- return std::make_shared<ga_instruction_reduction_opt0_2_unrolled<5>>
+ return
std::make_shared<ga_instruction_contraction_opt0_2_unrolled<5>>
(t, tc1, tc2, q2);
default:
- return std::make_shared<ga_instruction_reduction_opt0_2>
+ return std::make_shared<ga_instruction_contraction_opt0_2>
(t,tc1,tc2,n2,q2);
}
}
@@ -2776,66 +2996,66 @@ namespace getfem {
switch (q1) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<1,2>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<1,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<1,3>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<1,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<1,4>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<1,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt2_0_unrolled<1>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<1>>
(t, tc1, tc2, q1);
}
case 2:
switch (q1) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<2,2>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<2,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<2,3>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<2,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<2,4>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<2,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt2_0_unrolled<2>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<2>>
(t, tc1, tc2, q1);
}
case 3:
switch (q1) {
case 2:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<3,2>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<3,2>>
(t, tc1, tc2);
case 3:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<3,3>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<3,3>>
(t, tc1, tc2);
case 4:
return
- std::make_shared<ga_instruction_reduction_opt2_0_dunrolled<3,4>>
+
std::make_shared<ga_instruction_contraction_opt2_0_dunrolled<3,4>>
(t, tc1, tc2);
default :
- return
std::make_shared<ga_instruction_reduction_opt2_0_unrolled<3>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<3>>
(t, tc1, tc2, q1);
}
- return std::make_shared<ga_instruction_reduction_opt2_0_unrolled<3>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<3>>
(t, tc1, tc2, q1);
case 4:
- return std::make_shared<ga_instruction_reduction_opt2_0_unrolled<4>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<4>>
(t, tc1, tc2, q1);
case 5:
- return std::make_shared<ga_instruction_reduction_opt2_0_unrolled<5>>
+ return
std::make_shared<ga_instruction_contraction_opt2_0_unrolled<5>>
(t, tc1, tc2, q1);
default:
- return std::make_shared<ga_instruction_reduction_opt2_0>
+ return std::make_shared<ga_instruction_contraction_opt2_0>
(t,tc1,tc2, n1, q1);
}
}
@@ -2849,82 +3069,82 @@ namespace getfem {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,1>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,1>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,1>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 2 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,2>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,2>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,2>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 3 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,3>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,3>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,3>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 4 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,4>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,4>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,4>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 5 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,5>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,5>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,5>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 6 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,6>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,6>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,6>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 7 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,7>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,7>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,7>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 8 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,8>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,8>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,8>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 9 :
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,9>>(t, tc1, tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,9>>(t, tc1, tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,9>>(t, tc1, tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
case 10:
switch(n) {
case 2: return std::make_shared<ga_ins_red_d_unrolled<2,10>>(t, tc1,
tc2);
case 3: return std::make_shared<ga_ins_red_d_unrolled<3,10>>(t, tc1,
tc2);
case 4: return std::make_shared<ga_ins_red_d_unrolled<4,10>>(t, tc1,
tc2);
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return
ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
- default: return ga_instruction_reduction_switch(t_,tc1_,tc2_,n,to_clear);
+ default: return ga_instruction_contraction_switch(t_,tc1_,tc2_,n,to_clear);
}
}
// Performs Amij Bnj -> Cmni. To be optimized.
- struct ga_instruction_spec_reduction : public ga_instruction {
+ struct ga_instruction_spec_contraction : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type nn;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: specific reduction operation of "
+ GA_DEBUG_INFO("Instruction: specific contraction operation of "
"size " << nn);
size_type s1 = tc1.sizes()[0], s11 = tc1.size() / (s1*nn), s111 = s1*s11;
size_type s2 = tc2.sizes()[0];
@@ -2939,17 +3159,17 @@ namespace getfem {
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
- ga_instruction_spec_reduction(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_spec_contraction(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_) {}
};
// Performs Amik Bnjk -> Cmnij. To be optimized.
- struct ga_instruction_spec2_reduction : public ga_instruction {
+ struct ga_instruction_spec2_contraction : public ga_instruction {
base_tensor &t, &tc1, &tc2;
size_type nn;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: second specific reduction operation of "
+ GA_DEBUG_INFO("Instruction: second specific contraction operation of "
"size " << nn);
size_type s1 = tc1.sizes()[0], s11 = tc1.size() / (s1*nn), s111 = s1*s11;
size_type s2 = tc2.sizes()[0], s22 = tc2.size() / (s2*nn), s222 = s2*s22;
@@ -2965,7 +3185,7 @@ namespace getfem {
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
- ga_instruction_spec2_reduction(base_tensor &t_, base_tensor &tc1_,
+ ga_instruction_spec2_contraction(base_tensor &t_, base_tensor &tc1_,
base_tensor &tc2_, size_type n_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_) {}
};
@@ -4381,7 +4601,8 @@ namespace getfem {
pnode->node_type == GA_NODE_SPEC_FUNC ||
pnode->node_type == GA_NODE_CONSTANT ||
pnode->node_type == GA_NODE_ALLINDICES ||
- pnode->node_type == GA_NODE_RESHAPE) return;
+ pnode->node_type == GA_NODE_RESHAPE ||
+ pnode->node_type == GA_NODE_CONTRACT) return;
// cout << "compiling "; ga_print_node(pnode, cout); cout << endl;
@@ -4554,7 +4775,7 @@ namespace getfem {
case GA_NODE_PREDEF_FUNC: case GA_NODE_OPERATOR: case GA_NODE_SPEC_FUNC:
case GA_NODE_CONSTANT: case GA_NODE_ALLINDICES: case GA_NODE_ZERO:
- case GA_NODE_RESHAPE: case GA_NODE_INTERPOLATE_FILTER:
+ case GA_NODE_RESHAPE: case GA_NODE_CONTRACT: case
GA_NODE_INTERPOLATE_FILTER:
break;
case GA_NODE_X:
@@ -5310,16 +5531,17 @@ namespace getfem {
case GA_DOT: case GA_COLON: case GA_MULT:
{
- size_type tps1 = child0->tensor_proper_size();
- size_type tps2 = child1->tensor_proper_size();
- size_type s1 = (tps1 * tps2) / pnode->tensor_proper_size();
+ size_type tps0 = child0->tensor_proper_size();
+ size_type tps1 = child1->tensor_proper_size();
+ size_type s1 = (tps0 * tps1) / pnode->tensor_proper_size();
size_type s2 = size_type(round(sqrt(scalar_type(s1))));
pgai = pga_instruction();
- if (pnode->op_type == GA_DOT || pnode->op_type == GA_COLON ||
+ if ((pnode->op_type == GA_DOT && dim1 <= 1) ||
+ pnode->op_type == GA_COLON ||
(pnode->op_type == GA_MULT && dim0 == 4) ||
(pnode->op_type == GA_MULT && dim1 <= 1) ||
- child0->tensor().size() == 1 || child1->tensor().size() == 1) {
+ child0->tensor().size() == 1 || tps1 == 1) {
if (child0->tensor().size() == 1 && child1->tensor().size() == 1)
{
pgai = std::make_shared<ga_instruction_scalar_scalar_mult>
@@ -5336,7 +5558,7 @@ namespace getfem {
(pnode->tensor(), child0->tensor(), child1->tensor()[0]);
}
else if (pnode->test_function_type < 3) {
- if (child0->tensor_proper_size() == 1) {
+ if (tps0 == 1) {
if (is_uniform) // Unrolled instruction
pgai = ga_uniform_instruction_simple_tmult
(pnode->tensor(), child0->tensor(), child1->tensor());
@@ -5352,10 +5574,10 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_simple_tmult>
(pnode->tensor(), child1->tensor(), child0->tensor());
} else if (is_uniform) // Unrolled instruction
- pgai = ga_uniform_instruction_reduction_switch
+ pgai = ga_uniform_instruction_contraction_switch
(pnode->t, child0->t, child1->t, s2, tensor_to_clear);
else // Unrolled instruction
- pgai = ga_instruction_reduction_switch
+ pgai = ga_instruction_contraction_switch
(pnode->t, child0->t, child1->t, s2, tensor_to_clear);
}
} else {
@@ -5363,7 +5585,7 @@ namespace getfem {
child1->test_function_type == 3) {
if (child1->test_function_type == 3 ||
child1->tensor_proper_size() <= s2) {
- if (tps1 == 1) {
+ if (tps0 == 1) {
if (is_uniform) { // Unrolled instruction
pgai = ga_uniform_instruction_simple_tmult
(pnode->tensor(), child1->tensor(), child0->tensor());
@@ -5371,18 +5593,18 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_simple_tmult>
(pnode->tensor(), child1->tensor(), child0->tensor());
} else if (is_uniform) // Unrolled instruction
- pgai = ga_uniform_instruction_reduction_switch
+ pgai = ga_uniform_instruction_contraction_switch
(pnode->t, child0->t, child1->t, s2, tensor_to_clear);
else // Unrolled instruction
- pgai = ga_instruction_reduction_switch
+ pgai = ga_instruction_contraction_switch
(pnode->t, child0->t, child1->t, s2, tensor_to_clear);
} else
- pgai = std::make_shared<ga_instruction_spec_reduction>
+ pgai = std::make_shared<ga_instruction_spec_contraction>
(pnode->tensor(), child1->tensor(), child0->tensor(), s2);
} else if (child1->test_function_type == 0 ||
(child0->tensor_proper_size() == s2 &&
child1->tensor_proper_size() == s2)) {
- if (tps1 == 1) {
+ if (tps0 == 1) {
if (is_uniform) { // Unrolled instruction
pgai = ga_uniform_instruction_simple_tmult
(pnode->tensor(), child0->tensor(), child1->tensor());
@@ -5391,37 +5613,28 @@ namespace getfem {
(pnode->tensor(), child0->tensor(), child1->tensor());
} else {
if (is_uniform) // Unrolled instruction
- pgai = ga_uniform_instruction_reduction_switch
+ pgai = ga_uniform_instruction_contraction_switch
(pnode->t, child1->t, child0->t, s2, tensor_to_clear);
else // Unrolled instruction
- pgai = ga_instruction_reduction_switch
+ pgai = ga_instruction_contraction_switch
(pnode->t, child1->t, child0->t, s2, tensor_to_clear);
}
} else {
if (child0->tensor_proper_size() == s2)
- pgai = ga_uniform_instruction_reduction_switch
+ pgai = ga_uniform_instruction_contraction_switch
(pnode->t, child1->t, child0->t, s2, tensor_to_clear);
else if (child1->tensor_proper_size() == s2)
- pgai = std::make_shared<ga_instruction_spec_reduction>
+ pgai = std::make_shared<ga_instruction_spec_contraction>
(pnode->tensor(), child0->tensor(), child1->tensor(), s2);
else
- pgai = std::make_shared<ga_instruction_spec2_reduction>
+ pgai = std::make_shared<ga_instruction_spec2_contraction>
(pnode->tensor(), child0->tensor(), child1->tensor(), s2);
}
}
-
-
- } else { // GA_MULT
-
+ } else { // GA_MULT or GA_DOT for dim1 > 1
+ // and child1->tensor_proper_size() > 1
if (pnode->test_function_type < 3) {
- if (child1->tensor_proper_size() == 1) {
- if (is_uniform) // Unrolled instruction
- pgai = ga_uniform_instruction_simple_tmult
- (pnode->tensor(), child1->tensor(), child0->tensor());
- else
- pgai = std::make_shared<ga_instruction_simple_tmult>
- (pnode->tensor(), child1->tensor(), child0->tensor());
- } else if (child0->tensor_proper_size() == 1) {
+ if (tps0 == 1) {
if (is_uniform) // Unrolled instruction
pgai = ga_uniform_instruction_simple_tmult
(pnode->tensor(), child0->tensor(), child1->tensor());
@@ -5429,26 +5642,16 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_simple_tmult>
(pnode->tensor(), child0->tensor(), child1->tensor());
} else {
- if (dim0 == 2)
- pgai = std::make_shared<ga_instruction_matrix_mult>
- (pnode->tensor(), child0->tensor(), child1->tensor());
+ if (child1->test_function_type == 0)
+ pgai = std::make_shared<ga_instruction_matrix_mult>
+ (pnode->tensor(), child0->tensor(), child1->tensor(), s2);
+ else
+ pgai = std::make_shared<ga_instruction_matrix_mult_spec>
+ (pnode->tensor(), child0->tensor(), child1->tensor(),
+ s2, tps0/s2, tps1/s2);
}
} else {
- if (child1->tensor_proper_size() == 1) {
- if (child1->test_function_type == 0 ||
- child1->test_function_type == 1) {
- if (is_uniform) // Unrolled instruction
- pgai = ga_uniform_instruction_simple_tmult
- (pnode->tensor(), child1->tensor(), child0->tensor());
- else
- pgai = std::make_shared<ga_instruction_simple_tmult>
- (pnode->tensor(), child1->tensor(), child0->tensor());
- } else
- pgai = std::make_shared<ga_instruction_spec_tmult>
- (pnode->tensor(), child0->tensor(), child1->tensor(),
- child0->tensor_proper_size(),
- child1->tensor_proper_size());
- } else if (child0->tensor_proper_size() == 1) {
+ if (child0->tensor_proper_size() == 1) {
if (child0->test_function_type == 0 ||
child0->test_function_type == 1) {
if (is_uniform) // Unrolled instruction
@@ -5460,19 +5663,22 @@ namespace getfem {
} else
pgai = std::make_shared<ga_instruction_spec_tmult>
(pnode->tensor(), child1->tensor(), child0->tensor(),
- child1->tensor_proper_size(),
- child0->tensor_proper_size());
- } else if (dim0 == 2) {
- if (child1->test_function_type != 2)
+ tps1, tps0);
+ } else {
+ if (child1->test_function_type == 0)
pgai = std::make_shared<ga_instruction_matrix_mult>
- (pnode->tensor(), child0->tensor(), child1->tensor());
- else
- pgai = std::make_shared<ga_instruction_matrix_mult_spec>
- (pnode->tensor(), child0->tensor(), child1->tensor());
+ (pnode->tensor(), child0->tensor(), child1->tensor(), s2);
+ else if (child1->test_function_type == 2)
+ pgai = std::make_shared<ga_instruction_matrix_mult_spec>
+ (pnode->tensor(), child0->tensor(), child1->tensor(),
+ s2, tps0/s2, tps1/s2);
+ else
+ pgai = std::make_shared<ga_instruction_matrix_mult_spec2>
+ (pnode->tensor(), child0->tensor(), child1->tensor(),
+ s2, tps0/s2, tps1/s2);
}
}
}
- GMM_ASSERT1(pgai.get(), "Internal error");
rmi.instructions.push_back(std::move(pgai));
}
break;
@@ -5686,6 +5892,88 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_copy_tensor>(pnode->tensor(),
child1->tensor());
rmi.instructions.push_back(std::move(pgai));
+ } else if (child0->node_type == GA_NODE_CONTRACT) {
+ std::vector<size_type> ind(2), indsize(2);
+ pga_tree_node child2(0);
+ if (pnode->children.size() == 4)
+ { ind[0] = 2; ind[1] = 3; }
+ else if (pnode->children.size() == 5)
+ { ind[0] = 2; ind[1] = 4; child2 = pnode->children[3]; }
+ else if (pnode->children.size() == 7) {
+ ind.resize(4); indsize.resize(4);
+ ind[0] = 2; ind[1] = 3; ind[2] = 5; ind[3] = 6;
+ child2 = pnode->children[4];
+ }
+ size_type kk = 0, ll = 1;
+ for (size_type i = 1; i < pnode->children.size(); ++i) {
+ if (i == ind[kk]) {
+ ind[kk] = size_type(round(pnode->children[i]->tensor()[0]));
+ indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
+ ++kk;
+ } else ll = i;
+ }
+
+ if (pnode->children.size() == 4) {
+ size_type i1 = ind[0], i2 = ind[1];
+ if (i1 > i2) std::swap(i1, i2);
+ size_type ii2 = 1, ii3 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
+ if (i > i2) ii3 *= child1->tensor_proper_size(i);
+ }
+ pgai = std::make_shared<ga_instruction_contract_1_1>
+ (pnode->tensor(), child1->tensor(), indsize[0], ii2, ii3);
+ }
+ else if (pnode->children.size() == 5) {
+ // Particular cases should be detected (ii2=ii3=1 in particular).
+ size_type i1 = ind[0], i2 = ind[1];
+ size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1) ii2 *= child1->tensor_proper_size(i);
+ }
+ for (size_type i = 0; i < child2->tensor_order(); ++i) {
+ if (i < i2) ii3 *= child2->tensor_proper_size(i);
+ if (i > i2) ii4 *= child2->tensor_proper_size(i);
+ }
+ if (child1->test_function_type==1 && child2->test_function_type==2)
+ pgai = std::make_shared<ga_instruction_contract_2_1_rev>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ indsize[0], ii1, ii2, ii3, ii4);
+ else
+ pgai = std::make_shared<ga_instruction_contract_2_1>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ indsize[0], ii1, ii2, ii3, ii4);
+ }
+ else if (pnode->children.size() == 7) {
+ // Particular cases should be detected (ii2=ii3=1 in particular).
+ GMM_ASSERT1(false, "to be done !");
+ size_type i1 = ind[0], i2 = ind[1], i3 = ind[2], i4 = ind[3];
+ size_type nn1 = indsize[0], nn2 = indsize[1];
+ size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1, ii5 = 1, ii6 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
+ if (i > i2) ii3 *= child1->tensor_proper_size(i);
+ }
+ for (size_type i = 0; i < child2->tensor_order(); ++i) {
+ if (i < i3 && i < i4) ii4 *= child2->tensor_proper_size(i);
+ if ((i > i3 && i < i4) || (i > i4 && i < i3))
+ ii5 *= child2->tensor_proper_size(i);
+ if (i > i3 && i > i4) ii6 *= child2->tensor_proper_size(i);
+ }
+ if (i1 > i2)
+ { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
+ if (child1->test_function_type==1 && child2->test_function_type==2)
+ pgai = std::make_shared<ga_instruction_contract_2_2_rev>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6, i4 < i3);
+ else
+ pgai = std::make_shared<ga_instruction_contract_2_2>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6, i4 < i3);
+ }
+ rmi.instructions.push_back(std::move(pgai));
} else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
std::string name = child0->name;
diff --git a/src/getfem_generic_assembly_functions_and_operators.cc
b/src/getfem_generic_assembly_functions_and_operators.cc
index cb2c7e2..8f69e30 100644
--- a/src/getfem_generic_assembly_functions_and_operators.cc
+++ b/src/getfem_generic_assembly_functions_and_operators.cc
@@ -561,7 +561,6 @@ namespace getfem {
SPEC_OP.insert("element_K");
SPEC_OP.insert("element_B");
SPEC_OP.insert("Normal");
- SPEC_OP.insert("Reshape");
SPEC_OP.insert("Sym");
SPEC_OP.insert("Skew");
SPEC_OP.insert("Def");
@@ -573,6 +572,8 @@ namespace getfem {
SPEC_OP.insert("Xfem_plus");
SPEC_OP.insert("Xfem_minus");
SPEC_OP.insert("Print");
+ SPEC_OP.insert("Reshape");
+ SPEC_OP.insert("Contract");
SPEC_OP.insert("Diff");
}
diff --git a/src/getfem_generic_assembly_interpolation.cc
b/src/getfem_generic_assembly_interpolation.cc
index 4a088a9..fcff8f7 100644
--- a/src/getfem_generic_assembly_interpolation.cc
+++ b/src/getfem_generic_assembly_interpolation.cc
@@ -750,7 +750,7 @@ namespace getfem {
void init(const ga_workspace &/* workspace */) const {}
void finalize() const {}
- std::string expression(void) const { return "Id(meshdim)"; }
+ std::string expression(void) const { return "X"; }
int transform(const ga_workspace &/*workspace*/, const mesh &m_x,
fem_interpolation_context &ctx_x,
@@ -814,7 +814,7 @@ namespace getfem {
const std::string &/* interpolate_name */) const {}
void init(const ga_workspace &/* workspace */) const {}
void finalize() const {}
- std::string expression(void) const { return "Id(meshdim)"; }
+ std::string expression(void) const { return "X"; }
int transform(const ga_workspace &/*workspace*/, const mesh &m_x,
fem_interpolation_context &ctx_x,
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index d5a0731..6bd1c8d 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -294,7 +294,7 @@ namespace getfem {
case GA_NODE_PREDEF_FUNC: case GA_NODE_OPERATOR: case GA_NODE_SPEC_FUNC :
case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_ELT_SIZE:
case GA_NODE_ELT_K: case GA_NODE_ELT_B:
- case GA_NODE_NORMAL: case GA_NODE_RESHAPE:
+ case GA_NODE_NORMAL: case GA_NODE_RESHAPE: case GA_NODE_CONTRACT:
case GA_NODE_INTERPOLATE_X: case GA_NODE_INTERPOLATE_NORMAL:
pnode->test_function_type = 0; break;
@@ -1008,40 +1008,38 @@ namespace getfem {
break;
case GA_DOT:
- if (dim1 > 1)
- ga_throw_error(pnode->expr, pnode->pos, "The second argument of the "
- "dot product has to be a vector.")
- else {
+ {
size_type s0 = dim0 == 0 ? 1 : size0.back();
- size_type s1 = dim1 == 0 ? 1 : size1.back();
+ size_type s1 = dim1 == 0 ? 1 : size1[size1.size()-dim1];
if (s0 != s1) ga_throw_error(pnode->expr, pnode->pos, "Dot product "
"of expressions of different sizes ("
<< s0 << " != " << s1 << ").");
- if (child0->tensor_order() <= 1) pnode->symmetric_op = true;
+ if (dim0 <= 1 && dim1 <= 1) pnode->symmetric_op = true;
pnode->mult_test(child0, child1);
- if (dim0 > 1) {
+ if (dim0 > 1 || dim1 > 1) {
mi = pnode->t.sizes();
for (size_type i = 1; i < dim0; ++i)
mi.push_back(child0->tensor_proper_size(i-1));
+ for (size_type i = 1; i < dim1; ++i)
+ mi.push_back(child1->tensor_proper_size(i));
pnode->t.adjust_sizes(mi);
}
- if (all_cte) {
- pnode->node_type = GA_NODE_CONSTANT;
+ if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
+ gmm::clear(pnode->tensor().as_vector());
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ } else if (all_cte) {
gmm::clear(pnode->tensor().as_vector());
- size_type k = 0;
- for (size_type i = 0, j = 0; i < child0->tensor().size(); ++i) {
- pnode->tensor()[j] += child0->tensor()[i] * child1->tensor()[k];
- ++j; if (j == pnode->tensor().size()) { j = 0; ++k; }
- }
- GMM_ASSERT1(k == child1->tensor().size(), "Internal error");
+ size_type m0 = child0->tensor().size() / s0;
+ size_type m1 = child1->tensor().size() / s1;
+ for (size_type i = 0; i < m0; ++i)
+ for (size_type j = 0; j < m1; ++j)
+ for (size_type k = 0; k < s0; ++k)
+ pnode->tensor()[i*m1+j]
+ += child0->tensor()[i*s0+k] * child1->tensor()[k*m1+j];
+ pnode->node_type = GA_NODE_CONSTANT;
tree.clear_children(pnode);
- } else {
- if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
- gmm::clear(pnode->tensor().as_vector());
- pnode->node_type = GA_NODE_ZERO;
- tree.clear_children(pnode);
- }
}
}
break;
@@ -1070,7 +1068,12 @@ namespace getfem {
mi.push_back(child0->tensor_proper_size(i-2));
pnode->t.adjust_sizes(mi);
}
- if (all_cte) {
+
+ if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
+ gmm::clear(pnode->tensor().as_vector());
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ } else if (all_cte) {
pnode->node_type = GA_NODE_CONSTANT;
gmm::clear(pnode->tensor().as_vector());
size_type k = 0;
@@ -1080,12 +1083,6 @@ namespace getfem {
}
GMM_ASSERT1(k == child1->tensor().size(), "Internal error");
tree.clear_children(pnode);
- } else {
- if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
- gmm::clear(pnode->tensor().as_vector());
- pnode->node_type = GA_NODE_ZERO;
- tree.clear_children(pnode);
- }
}
}
break;
@@ -1448,7 +1445,12 @@ namespace getfem {
}
if (!(name.compare("Reshape"))) {
pnode->node_type = GA_NODE_RESHAPE;
- pnode->init_vector_tensor(meshdim);
+ pnode->init_scalar_tensor(0);
+ break;
+ }
+ if (!(name.compare("Contract"))) {
+ pnode->node_type = GA_NODE_CONTRACT;
+ pnode->init_scalar_tensor(0);
break;
}
@@ -1674,8 +1676,10 @@ namespace getfem {
break;
case GA_NODE_PARAMS:
+
+ // Diff operator
if (child0->node_type == GA_NODE_NAME) {
- if (child0->name.compare("Diff") == 0) { // Diff operator
+ if (child0->name.compare("Diff") == 0) {
if (pnode->children.size() != 3)
ga_throw_error(pnode->expr, child0->pos,
@@ -1708,7 +1712,10 @@ namespace getfem {
pnode = child1;
} else
ga_throw_error(pnode->expr, child0->pos, "Unknown special operator");
- } else if (child0->node_type == GA_NODE_X) {
+ }
+
+ // X (current coordinates on the mesh)
+ else if (child0->node_type == GA_NODE_X) {
child0->init_scalar_tensor(0);
if (pnode->children.size() != 2)
ga_throw_error(pnode->expr, child1->pos,
@@ -1725,8 +1732,10 @@ namespace getfem {
<< meshdim);
tree.replace_node_by_child(pnode, 0);
pnode = child0;
+ }
- } else if (child0->node_type == GA_NODE_RESHAPE) {
+ // Reshape operation
+ else if (child0->node_type == GA_NODE_RESHAPE) {
if (pnode->children.size() < 3)
ga_throw_error(pnode->expr, child1->pos,
"Not enough parameters for Reshape");
@@ -1768,9 +1777,199 @@ namespace getfem {
pnode->node_type = GA_NODE_ZERO;
tree.clear_children(pnode);
}
- } else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
+ }
+
+ // Tensor contraction operator
+ else if (child0->node_type == GA_NODE_CONTRACT) {
+ std::vector<size_type> ind(2), indsize(2);
+ if (pnode->children.size() == 4)
+ { ind[0] = 2; ind[1] = 3; }
+ else if (pnode->children.size() == 5)
+ { ind[0] = 2; ind[1] = 4; }
+ else if (pnode->children.size() == 7) {
+ ind.resize(4); indsize.resize(4);
+ ind[0] = 2; ind[1] = 3; ind[2] = 5; ind[3] = 6;
+ }
+
+ size_type order = 0, kk = 0, ll = 1; // Indices extraction and controls
+ for (size_type i = 1; i < pnode->children.size(); ++i) {
+ if (i == ind[kk]) {
+ if (pnode->children[i]->node_type != GA_NODE_CONSTANT ||
+ pnode->children[i]->tensor().size() != 1)
+ ga_throw_error(pnode->children[i]->expr, pnode->children[i]->pos,
+ "Invalid parameter for Contract. "
+ "Should be an index number.");
+ ind[kk] = size_type(round(pnode->children[i]->tensor()[0]));
+ indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
+ order = pnode->children[ll]->tensor_order();
+ if (ind[kk] < 1 || ind[kk] > order)
+ ga_throw_error(pnode->children[i]->expr, pnode->children[i]->pos,
+ "Parameter out of range for Contract (should be "
+ "between 1 and " << order << ")");
+ if (kk >= ind.size()/2 && indsize[kk] != indsize[kk-ind.size()/2])
+ ga_throw_error(child0->expr, child0->pos,
+ "Invalid parameters for Contract. Cannot "
+ "contract on indices of different sizes");
+ ++kk;
+ } else ll = i;
+ }
+
+ if (pnode->children.size() == 4) {
+ // Contraction of a single tensor on a single pair of indices
+ pnode->test_function_type = child1->test_function_type;
+ pnode->name_test1 = child1->name_test1;
+ pnode->name_test2 = child1->name_test2;
+ pnode->interpolate_name_test1 = child1->interpolate_name_test1;
+ pnode->interpolate_name_test2 = child1->interpolate_name_test2;
+ pnode->qdim1 = child1->qdim1;
+ pnode->qdim2 = child1->qdim2;
+
+ size_type i1 = ind[0], i2 = ind[1];
+ if (i1 == i2)
+ ga_throw_error(child0->expr, child0->pos,
+ "Invalid parameters for Contract. Repeated index");
+
+ mi.resize(0);
+ for (size_type i = 0; i < pnode->nb_test_functions(); ++i)
+ mi.push_back(size1[i]);
+ for (size_type i = 0; i < order; ++i)
+ if (i != i1 && i != i2)
+ mi.push_back(child1->tensor_proper_size(i));
+ pnode->t.adjust_sizes(mi);
+
+ if (child1->node_type == GA_NODE_CONSTANT) {
+
+ if (i1 > i2) std::swap(i1, i2);
+ size_type ii1 = 1, ii2 = 1, ii3 = 1;
+ size_type nn = indsize[0];
+ for (size_type i = 0; i < order; ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
+ if (i > i2) ii3 *= child1->tensor_proper_size(i);
+ }
+
+ auto it = pnode->tensor().begin();
+ for (size_type i = 0; i < ii3; ++i)
+ for (size_type j = 0; j < ii2; ++j)
+ for (size_type k = 0; k < ii1; ++k, ++it) {
+ *it = scalar_type(0);
+ size_type pre_ind = k+j*ii1*nn+i*ii1*nn*ii2*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += child1->tensor()[pre_ind+n*ii1+n*ii1*nn*ii2];
+ }
+
+ pnode->node_type = GA_NODE_CONSTANT;
+ tree.clear_children(pnode);
+ } else if (child1->node_type == GA_NODE_ZERO) {
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ }
+
+ } else if (pnode->children.size() == 5) {
+ // Contraction of two tensors on a single pair of indices
+ pga_tree_node child2 = pnode->children[3];
+ pnode->mult_test(child1, child2);
+
+ size_type i1 = ind[0], i2 = ind[1];
+ mi = pnode->tensor().sizes();
+ for (size_type i = 0; i < dim1; ++i)
+ if (i != i1) mi.push_back(child1->tensor_proper_size(i));
+ for (size_type i = 0; i < order; ++i)
+ if (i != i2) mi.push_back(child2->tensor_proper_size(i));
+ pnode->t.adjust_sizes(mi);
+
+ if (child1->node_type == GA_NODE_CONSTANT &&
+ child2->node_type == GA_NODE_CONSTANT) {
+ size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1;
+ size_type nn = indsize[0];
+ for (size_type i = 0; i < dim1; ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1) ii2 *= child1->tensor_proper_size(i);
+ }
+ for (size_type i = 0; i < order; ++i) {
+ if (i < i2) ii3 *= child2->tensor_proper_size(i);
+ if (i > i2) ii4 *= child2->tensor_proper_size(i);
+ }
+
+ auto it = pnode->tensor().begin();
+ for (size_type i = 0; i < ii4; ++i)
+ for (size_type j = 0; j < ii3; ++j)
+ for (size_type k = 0; k < ii2; ++k)
+ for (size_type l = 0; l < ii1; ++l, ++it) {
+ *it = scalar_type(0);
+ for (size_type n = 0; n < nn; ++n)
+ *it += child1->tensor()[l+n*ii1+k*ii1*nn]
+ * child2->tensor()[j+n*ii3+i*ii3*nn];
+ }
+
+ } else if (child1->tensor_is_zero() || child2->tensor_is_zero()) {
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ }
+
+ } else if (pnode->children.size() == 7) {
+ // Contraction of two tensors on two pairs of indices
+ pga_tree_node child2 = pnode->children[4];
+ pnode->mult_test(child1, child2);
+ if (ind[0] == ind[1] || ind[2] == ind[3])
+ ga_throw_error(child0->expr, child0->pos,
+ "Invalid parameters for Contract. Repeated index.");
+
+ size_type i1 = ind[0], i2 = ind[1], i3 = ind[2], i4 = ind[3];
+ mi = pnode->tensor().sizes();
+ for (size_type i = 0; i < dim1; ++i)
+ if (i != i1 && i != i2) mi.push_back(child1->tensor_proper_size(i));
+ for (size_type i = 0; i < order; ++i)
+ if (i != i3 && i != i4) mi.push_back(child2->tensor_proper_size(i));
+ pnode->t.adjust_sizes(mi);
+
+ if (child1->tensor_is_zero() || child2->tensor_is_zero()) {
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ } else if (child1->node_type == GA_NODE_CONSTANT &&
+ child2->node_type == GA_NODE_CONSTANT) {
+ size_type nn1 = indsize[0], nn2 = indsize[1];
+ if (i1 > i2)
+ { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
+
+ size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1, ii5 = 1, ii6 = 1;
+ for (size_type i = 0; i < dim1; ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
+ if (i > i2) ii3 *= child1->tensor_proper_size(i);
+ }
+ for (size_type i = 0; i < order; ++i) {
+ if (i < i3 && i < i4) ii4 *= child2->tensor_proper_size(i);
+ if ((i > i3 && i < i4) || (i > i4 && i < i3))
+ ii5 *= child2->tensor_proper_size(i);
+ if (i > i3 && i > i4) ii6 *= child2->tensor_proper_size(i);
+ }
+
+ auto it = pnode->tensor().begin();
+ size_type m1 = ii4, m2 = ii4*nn1*ii5;
+ if (i3 < i4) std::swap(m1, m2);
+ for (size_type i = 0; i < ii6; ++i)
+ for (size_type j = 0; j < ii5; ++j)
+ for (size_type k = 0; k < ii4; ++k)
+ for (size_type l = 0; l < ii3; ++l)
+ for (size_type m = 0; m < ii2; ++m)
+ for (size_type p = 0; p < ii1; ++p, ++it) {
+ *it = scalar_type(0);
+ size_type q1 = p+m*ii1*nn1+l*ii1*nn1*ii2*nn2;
+ size_type q2 = k+j*ii4*nn1+i*ii4*nn1*ii5*nn2;
+ for (size_type n1 = 0; n1 < nn1; ++n1)
+ for (size_type n2 = 0; n2 < nn2; ++n2)
+ *it += child1->tensor()[q1+n1*ii1+n2*ii1*nn1*ii2]
+ * child2->tensor()[q2+n1*m1+n2*m2];
+ }
+ }
+ } else
+ ga_throw_error(pnode->expr, child1->pos,
+ "Wrong number of parameters for Contract");
+ }
- // Evaluation of a predefined function
+ // Evaluation of a predefined function
+ else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
for (size_type i = 1; i < pnode->children.size(); ++i)
ga_valid_operand(pnode->children[i]);
@@ -1838,9 +2037,10 @@ namespace getfem {
}
tree.clear_children(pnode);
}
- } else if (child0->node_type == GA_NODE_SPEC_FUNC) {
+ }
- // Special constant functions: meshdim, qdim(u) ...
+ // Special constant functions: meshdim, qdim(u) ...
+ else if (child0->node_type == GA_NODE_SPEC_FUNC) {
for (size_type i = 1; i < pnode->children.size(); ++i)
ga_valid_operand(pnode->children[i]);
@@ -1888,14 +2088,16 @@ namespace getfem {
pnode->init_scalar_tensor(scalar_type(1));
else {
pnode->init_matrix_tensor(n,n);
+ gmm::clear(pnode->tensor().as_vector());
for (int i = 0; i < n; ++i) pnode->tensor()(i,i) = scalar_type(1);
}
} else ga_throw_error(pnode->expr, pnode->children[0]->pos,
"Unknown special function.");
tree.clear_children(pnode);
- } else if (child0->node_type == GA_NODE_OPERATOR) {
+ }
- // Call to a nonlinear operator
+ // Nonlinear operator call
+ else if (child0->node_type == GA_NODE_OPERATOR) {
for (size_type i = 1; i < pnode->children.size(); ++i)
ga_valid_operand(pnode->children[i]);
@@ -1962,13 +2164,11 @@ namespace getfem {
tree.clear_children(pnode);
}
}
- } else {
+ }
- // Access to components of a tensor
+ // Access to components of a tensor
+ else {
all_cte = (child0->node_type == GA_NODE_CONSTANT);
- // cout << "child0->tensor_order() = " << child0->tensor_order();
- // cout << endl << "child0->t.sizes() = "
- // << child0->t.sizes() << endl;
if (pnode->children.size() != child0->tensor_order() + 1)
ga_throw_error(pnode->expr, pnode->pos, "Bad number of indices.");
for (size_type i = 1; i < pnode->children.size(); ++i)
@@ -2006,7 +2206,11 @@ namespace getfem {
pnode->qdim1 = child0->qdim1;
pnode->qdim2 = child0->qdim2;
- if (all_cte) {
+ if (child0->tensor_is_zero()) {
+ gmm::clear(pnode->tensor().as_vector());
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ } else if (all_cte) {
pnode->node_type = GA_NODE_CONSTANT;
for (bgeot::multi_index mi3(mi2.size()); !mi3.finished(mi2);
mi3.incrementation(mi2)) {
@@ -2016,12 +2220,6 @@ namespace getfem {
pnode->tensor()(mi3) = pnode->children[0]->tensor()(mi1);
}
tree.clear_children(pnode);
- } else {
- if (child0->tensor_is_zero() || child1->tensor_is_zero()) {
- gmm::clear(pnode->tensor().as_vector());
- pnode->node_type = GA_NODE_ZERO;
- tree.clear_children(pnode);
- }
}
}
break;
@@ -2073,7 +2271,7 @@ namespace getfem {
void ga_extract_factor(ga_tree &result_tree, pga_tree_node pnode,
- pga_tree_node &new_pnode) {
+ pga_tree_node &new_pnode) {
result_tree.clear();
result_tree.copy_node(pnode, 0, result_tree.root);
new_pnode = result_tree.root;
@@ -2106,6 +2304,7 @@ namespace getfem {
result_tree.insert_node(result_tree.root, pnode->node_type);
result_tree.root->op_type = pnode->op_type;
result_tree.root->pos = pnode->pos;
+ result_tree.root->expr = pnode->expr;
break;
case GA_DOT: case GA_MULT: case GA_COLON: case GA_TMULT:
case GA_DOTMULT: case GA_DIV: case GA_DOTDIV:
@@ -2113,6 +2312,7 @@ namespace getfem {
result_tree.insert_node(result_tree.root, pnode->node_type);
result_tree.root->op_type = pnode->op_type;
result_tree.root->pos = pnode->pos;
+ result_tree.root->expr = pnode->expr;
result_tree.root->children.resize(2, nullptr);
if (child0 == pnode_child) {
result_tree.copy_node(child1, result_tree.root,
@@ -2129,25 +2329,42 @@ namespace getfem {
break;
case GA_NODE_PARAMS:
- GMM_ASSERT1(child0->node_type == GA_NODE_RESHAPE, "Cannot extract a "
- "factor which is a parameter of a nonlinear "
- "operator/function");
- GMM_ASSERT1(child1 == pnode_child, "Cannot extract a factor of a "
- "Reshape size parameter");
- // Copy of the term and other children
- result_tree.insert_node(result_tree.root, pnode->node_type);
- result_tree.root->pos = pnode->pos;
- result_tree.root->children.resize(pnode->children.size(), nullptr);
- std::swap(result_tree.root->children[1], result_tree.root->children[0]);
- for (size_type i = 0; i < pnode->children.size(); ++i)
- if (i != 1)
- result_tree.copy_node(pnode->children[i], result_tree.root,
- result_tree.root->children[i]);
+ if (child0->node_type == GA_NODE_RESHAPE) {
+ GMM_ASSERT1(child1 == pnode_child, "Cannot extract a factor of a "
+ "Reshape size parameter");
+ // Copy of the term and other children
+ result_tree.insert_node(result_tree.root, pnode->node_type);
+ result_tree.root->pos = pnode->pos;
+ result_tree.root->expr = pnode->expr;
+ result_tree.root->children.resize(pnode->children.size(), nullptr);
+
std::swap(result_tree.root->children[1],result_tree.root->children[0]);
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ if (i != 1)
+ result_tree.copy_node(pnode->children[i], result_tree.root,
+ result_tree.root->children[i]);
+ } else if (child0->node_type == GA_NODE_CONTRACT) {
+ // Copy of the term and other children
+ result_tree.insert_node(result_tree.root, pnode->node_type);
+ result_tree.root->pos = pnode->pos;
+ result_tree.root->expr = pnode->expr;
+ result_tree.root->children.resize(pnode->children.size(), nullptr);
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ if (pnode->children[i] == pnode_child)
+ std::swap(result_tree.root->children[i],
+ result_tree.root->children[0]);
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ if (pnode->children[i] != pnode_child)
+ result_tree.copy_node(pnode->children[i], result_tree.root,
+ result_tree.root->children[i]);
+ } else
+ GMM_ASSERT1(false, "Cannot extract a factor which is a parameter "
+ "of a nonlinear operator/function");
break;
case GA_NODE_C_MATRIX:
result_tree.insert_node(result_tree.root, pnode->node_type);
result_tree.root->pos = pnode->pos;
+ result_tree.root->expr = pnode->expr;
result_tree.root->children.resize(pnode->children.size());
for (size_type i = 0; i < pnode->children.size(); ++i)
if (pnode_child == pnode->children[i]) {
@@ -2157,8 +2374,8 @@ namespace getfem {
for (size_type i = 0; i < pnode->children.size(); ++i) {
if (pnode_child == pnode->children[i]) {
- pnode->children[i] = new ga_tree_node(GA_NODE_ZERO, pnode->pos,
- pnode->expr);
+ pnode->children[i]
+ = new ga_tree_node(GA_NODE_ZERO, pnode->pos, pnode->expr);
pnode->children[i]->init_scalar_tensor(scalar_type(0));
pnode->children[i]->parent = pnode;
}
@@ -2177,6 +2394,7 @@ namespace getfem {
result_tree.insert_node(result_tree.root, GA_NODE_OP);
result_tree.root->op_type = GA_UNARY_MINUS;
result_tree.root->pos = pnode->children[0]->pos;
+ result_tree.root->expr = pnode->children[0]->expr;
}
}
@@ -2203,8 +2421,8 @@ namespace getfem {
case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST: case GA_NODE_PREDEF_FUNC:
case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST: case GA_NODE_RESHAPE:
- case GA_NODE_ELT_SIZE: case GA_NODE_ELT_K: case GA_NODE_ELT_B:
- case GA_NODE_CONSTANT: case GA_NODE_X:
+ case GA_NODE_CONTRACT: case GA_NODE_ELT_SIZE: case GA_NODE_ELT_K:
+ case GA_NODE_ELT_B: case GA_NODE_CONSTANT: case GA_NODE_X:
case GA_NODE_NORMAL: case GA_NODE_OPERATOR:
is_constant = true; break;
@@ -2292,6 +2510,12 @@ namespace getfem {
case GA_NODE_PARAMS:
if (child0->node_type == GA_NODE_RESHAPE) {
is_constant = child_1_is_constant;
+ } else if (child0->node_type == GA_NODE_CONTRACT) {
+ for (size_type i = 1; i < pnode->children.size(); ++i) {
+ if (!(ga_node_extract_constant_term(tree, pnode->children[i],
+ workspace, m)))
+ { is_constant = false; break; }
+ }
} else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
for (size_type i = 1; i < pnode->children.size(); ++i) {
if (!(ga_node_extract_constant_term(tree, pnode->children[i],
@@ -2356,7 +2580,7 @@ namespace getfem {
// A*Div_Test_u --> A*Normal if A is a scalar
// Div_Test_u*A --> Normal*A if A is a scalar
// A*(Grad_Test_u)' --> (A)'*Normal if A is a matrix
- // intercaled scalar multplications and divisions are taken into account
+ // intercaled scalar multiplications and divisions are taken into account
while (nnew_pnode->parent) {
pga_tree_node previous_node = nnew_pnode;
nnew_pnode = nnew_pnode->parent;
@@ -2484,7 +2708,6 @@ namespace getfem {
ga_extract_Neumann_term(tree, varname, workspace,
pnode->children[i], result);
-
switch (pnode->node_type) {
case GA_NODE_DIVERG_TEST: case GA_NODE_GRAD_TEST:
case GA_NODE_ELEMENTARY_GRAD_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
@@ -2893,6 +3116,30 @@ namespace getfem {
if (child0->node_type == GA_NODE_RESHAPE) {
ga_node_derivation(tree, workspace, m, pnode->children[1],
varname, interpolatename, order);
+ } else if (child0->node_type == GA_NODE_CONTRACT) {
+
+ if (pnode->children.size() == 4) {
+ ga_node_derivation(tree, workspace, m, pnode->children[1],
+ varname, interpolatename, order);
+ } else if (pnode->children.size() == 5 || pnode->children.size() == 7) {
+ size_type n2 = (pnode->children.size()==5) ? 3 : 4;
+ pga_tree_node child2 = pnode->children[n2];
+
+ if (mark1 && child2->marked) {
+ tree.duplicate_with_addition(pnode);
+ ga_node_derivation(tree, workspace, m, child1, varname,
+ interpolatename, order);
+ ga_node_derivation(tree, workspace, m,
+ pnode->parent->children[1]->children[n2],
+ varname, interpolatename, order);
+ } else if (mark1) {
+ ga_node_derivation(tree, workspace, m, child1, varname,
+ interpolatename, order);
+ } else
+ ga_node_derivation(tree, workspace, m, child2, varname,
+ interpolatename, order);
+
+ } else GMM_ASSERT1(false, "internal error");
} else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
std::string name = child0->name;
ga_predef_function_tab::const_iterator it =
PREDEF_FUNCTIONS.find(name);
@@ -3066,7 +3313,7 @@ namespace getfem {
pnode2->children[0]->der1 = i;
tree.insert_node(pnode2, GA_NODE_OP);
pga_tree_node pnode_op = pnode2->parent;
- // calcul de l'ordre de reduction
+ // Computation of contraction order
size_type red = pnode->children[i]->tensor_order();
switch (red) {
case 0 : pnode_op->op_type = GA_MULT; break;
@@ -3200,12 +3447,12 @@ namespace getfem {
tree.duplicate_with_operation(pnode, GA_COLON);
child0 = pnode; pnode = pnode->parent; child1 = pnode->children[1];
child1->init_matrix_tensor(meshdim, meshdim);
+ gmm::clear(pnode->tensor().as_vector());
for (size_type i = 0; i < meshdim; ++i)
child1->tensor()(i,i) = scalar_type(1);
child1->node_type = GA_NODE_CONSTANT;
break;
-
case GA_NODE_INTERPOLATE_HESS_TEST:
case GA_NODE_INTERPOLATE_HESS:
GMM_ASSERT1(false, "Sorry, cannot derive a hessian once more");
@@ -3214,265 +3461,179 @@ namespace getfem {
case GA_NODE_INTERPOLATE_VAL:
case GA_NODE_INTERPOLATE_GRAD:
case GA_NODE_INTERPOLATE_DIVERG:
- {
- bool is_val(pnode->node_type == GA_NODE_INTERPOLATE_VAL);
- bool is_grad(pnode->node_type == GA_NODE_INTERPOLATE_GRAD);
- bool is_diverg(pnode->node_type == GA_NODE_INTERPOLATE_DIVERG);
-
+ case GA_NODE_INTERPOLATE_VAL_TEST:
+ case GA_NODE_INTERPOLATE_GRAD_TEST:
+ case GA_NODE_INTERPOLATE_DIVERG_TEST:
+ case GA_NODE_INTERPOLATE_X:
+ {
std::string &tname = pnode->interpolate_name;
- std::string expr_trans =
- workspace.interpolate_transformation(tname)->expression();
+ auto ptrans = workspace.interpolate_transformation(tname);
+ std::string expr_trans = ptrans->expression();
if (expr_trans.size() == 0)
GMM_ASSERT1(false, "Sorry, the gradient of tranformation "
<< tname << " cannot be calculated. "
"The gradient computation is available only for "
"transformations having an explicit expression");
+ ga_tree trans_tree;
+ ga_read_string(expr_trans, trans_tree, workspace.macro_dictionnary());
+ ga_semantic_analysis(trans_tree, workspace, m,
+ ref_elt_dim_of_mesh(m), false, false, 1);
+ if (trans_tree.root) {
+ ga_node_grad(trans_tree, workspace, m, trans_tree.root);
+ ga_semantic_analysis(trans_tree, workspace, m,
+ ref_elt_dim_of_mesh(m), false, false, 1);
- // bool ivar = (pnode->name.compare(varname) == 0 &&
- // pnode->interpolate_name.compare(interpolatename) == 0);
- // bool itrans = !ivar;
- // if (!itrans) {
- // std::set<var_trans_pair> vars;
- // workspace.interpolate_transformation(pnode->interpolate_name)
- // ->extract_variables(workspace, vars, true, m,
- // pnode->interpolate_name);
- // for (const var_trans_pair &var : vars) {
- // if (var.varname.compare(varname) == 0 &&
- // var.transname.compare(interpolatename) == 0)
- // itrans = true;
- // }
- // }
-
- pga_tree_node pnode_trans = pnode;
- tree.duplicate_with_operation(pnode_trans, GA_MULT);
- pga_tree_node pnode_grad_trans = pnode_trans->parent->children[1];
-
- if (is_val) pnode_trans->node_type = GA_NODE_INTERPOLATE_GRAD;
- if (is_grad) pnode_trans->node_type = GA_NODE_INTERPOLATE_HESS;
- // if (is_diverg) ??;
-
- // Pour le Hessien, il s'agit vraiment d'une multiplication ? ou d'un
"." ?
-
- // Calcul et insertion du gradient de la transformation
-
-// if (ivar) {
-// mi.resize(1); mi[0] = 2;
-// for (size_type i = 0; i < pnode->tensor_order(); ++i)
-// mi.push_back(pnode->tensor_proper_size(i));
-// pnode->t.adjust_sizes(mi);
-// if (is_val) // --> t(Qmult*ndof,Qmult*target_dim)
-// pnode->node_type = GA_NODE_INTERPOLATE_VAL_TEST;
-// else if (is_grad) // --> t(Qmult*ndof,Qmult*target_dim,N)
-// pnode->node_type = GA_NODE_INTERPOLATE_GRAD_TEST;
-// else if (is_hess) // --> t(Qmult*ndof,Qmult*target_dim,N,N)
-// pnode->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
-// else if (is_diverg) // --> t(Qmult*ndof)
-// pnode->node_type = GA_NODE_INTERPOLATE_DIVERG_TEST;
-// pnode->test_function_type = order;
-// }
-
-// if (itrans) {
-// const mesh_fem *mf = workspace.associated_mf(pnode_trans->name);
-// size_type q = workspace.qdim(pnode_trans->name);
-// size_type n = mf->linked_mesh().dim();
-// bgeot::multi_index mii = workspace.qdims(pnode_trans->name);
-
-// if (is_val) // --> t(target_dim*Qmult,N)
-// pnode_trans->node_type = GA_NODE_INTERPOLATE_GRAD;
-// else if (is_grad || is_diverg) // --> t(target_dim*Qmult,N,N)
-// pnode_trans->node_type = GA_NODE_INTERPOLATE_HESS;
-
-// if (n > 1) {
-// if (q == 1 && mii.size() <= 1) { mii.resize(1); mii[0] = n; }
-// else mii.push_back(n);
-
-// if (is_grad || is_diverg) mii.push_back(n);
-// }
-// pnode_trans->t.adjust_sizes(mii);
-// tree.duplicate_with_operation(pnode_trans,
-// (n > 1) ? GA_DOT : GA_MULT);
-// pga_tree_node pnode_der = pnode_trans->parent->children[1];
-// pnode_der->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
-// if (n == 1)
-// pnode_der->init_vector_tensor(2);
-// else
-// pnode_der->init_matrix_tensor(2, n);
-// pnode_der->test_function_type = order;
-// pnode_der->name = varname;
-// pnode_der->interpolate_name_der = pnode_der->interpolate_name;
-// pnode_der->interpolate_name = interpolatename;
-
-// if (is_diverg) { // --> t(Qmult*ndof)
-// tree.insert_node(pnode_trans->parent, GA_NODE_OP);
-// pga_tree_node pnode_tr = pnode_trans->parent->parent;
-// pnode_tr->op_type = GA_TRACE;
-// pnode_tr->init_vector_tensor(2);
-// // pnode_tr->test_function_type = order;
-// // pnode_tr->name_test1 = pnode_trans->name_test1;
-// // pnode_tr->name_test2 = pnode_trans->name_test2;
-// }
-// }
- }
- break;
- #ifdef continue_here
-
- case GA_NODE_INTERPOLATE_VAL_TEST:
- case GA_NODE_INTERPOLATE_GRAD_TEST:
- case GA_NODE_INTERPOLATE_DIVERG_TEST:
- {
- bool is_val(pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST);
- bool is_grad(pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST);
- bool is_diverg(pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST);
-
- pga_tree_node pnode_trans = pnode;
- const mesh_fem *mf = workspace.associated_mf(pnode_trans->name);
- size_type q = workspace.qdim(pnode_trans->name);
- size_type n = mf->linked_mesh().dim();
- bgeot::multi_index mii = workspace.qdims(pnode_trans->name);
- if (is_val) // --> t(Qmult*ndof,Qmult*target_dim,N)
- pnode_trans->node_type = GA_NODE_INTERPOLATE_GRAD_TEST;
- else if (is_grad || is_diverg) // -->
t(Qmult*ndof,Qmult*target_dim,N,N)
- pnode_trans->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
-
- if (q == 1 && mii.size() <= 1) { mii.resize(1); mii[0] = 2; }
- else mii.insert(mii.begin(), 2);
+ GMM_ASSERT1(trans_tree.root->tensor().sizes().size() == 2,
+ "Problem in transformation" << tname);
+ size_type trans_dim = trans_tree.root->tensor().sizes()[1];
- if (n > 1) {
- mii.push_back(n);
- if (is_grad || is_diverg) mii.push_back(n);
- }
- pnode_trans->t.adjust_sizes(mii);
- // pnode_trans->test_function_type = order;
- tree.duplicate_with_operation(pnode_trans,
- (n > 1 ? GA_DOT : GA_MULT));
- pga_tree_node pnode_der = pnode_trans->parent->children[1];
- pnode_der->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
- if (n == 1)
- pnode_der->init_vector_tensor(2);
- else
- pnode_der->init_matrix_tensor(2, n);
- pnode_der->test_function_type = order;
- pnode_der->name = varname;
- pnode_der->interpolate_name_der = pnode_der->interpolate_name;
- pnode_der->interpolate_name = interpolatename;
+ tree.duplicate_with_operation(pnode, GA_DOT);
+
+ if (pnode->node_type == GA_NODE_INTERPOLATE_VAL)
+ pnode->node_type = GA_NODE_INTERPOLATE_GRAD;
+ if (pnode->node_type == GA_NODE_INTERPOLATE_GRAD ||
+ pnode->node_type == GA_NODE_INTERPOLATE_DIVERG)
+ pnode->node_type = GA_NODE_INTERPOLATE_HESS;
+ if (pnode->node_type == GA_NODE_INTERPOLATE_VAL_TEST)
+ pnode->node_type = GA_NODE_INTERPOLATE_GRAD_TEST;
+ if (pnode->node_type == GA_NODE_INTERPOLATE_GRAD_TEST ||
+ pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST)
+ pnode->node_type = GA_NODE_INTERPOLATE_HESS_TEST;
+
+ if (pnode->node_type == GA_NODE_INTERPOLATE_DIVERG ||
+ pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST) {
+ tree.duplicate_with_operation(pnode, GA_COLON);
+ child0 = pnode; pnode = pnode->parent; child1 = pnode->children[1];
+ child1->init_matrix_tensor(trans_dim, trans_dim);
+ gmm::clear(pnode->tensor().as_vector());
+ for (size_type i = 0; i < trans_dim; ++i)
+ child1->tensor()(i,i) = scalar_type(1);
+ child1->node_type = GA_NODE_CONSTANT;
+ }
- if (is_diverg) { // --> t(Qmult*ndof)
- tree.insert_node(pnode_trans->parent, GA_NODE_OP);
- pga_tree_node pnode_tr = pnode_trans->parent->parent;
- pnode_tr->op_type = GA_TRACE;
- pnode_tr->init_vector_tensor(2);
-// pnode_tr->test_function_type = order;
-// pnode_tr->name_test1 = pnode_trans->name_test1;
-// pnode_tr->name_test2 = pnode_trans->name_test2;
- }
+ if (pnode->node_type == GA_NODE_INTERPOLATE_X) {
+ pnode->node_type = GA_NODE_CONSTANT;
+ pnode->init_matrix_tensor(trans_dim, trans_dim);
+ gmm::clear(pnode->tensor().as_vector());
+ for (size_type i = 0; i < trans_dim; ++i)
+ pnode->tensor()(i,i) = scalar_type(1);
+ }
+
+ tree.clear_node_rec(pnode->parent->children[1]);
+ pnode->parent->children[1] = nullptr;
+ tree.copy_node(trans_tree.root, pnode->parent,
+ pnode->parent->children[1]);
+ } else {
+ pnode->node_type = GA_NODE_ZERO;
+ gmm::clear(pnode->tensor().as_vector());
+ }
}
break;
-
-
- case GA_NODE_INTERPOLATE_X:
- {
- size_type n = m.dim();
- pga_tree_node pnode_der = pnode;
- pnode_der->node_type = GA_NODE_INTERPOLATE_DERIVATIVE;
- if (n == 1)
- pnode_der->init_vector_tensor(2);
- else
- pnode_der->init_matrix_tensor(2, n);
- pnode_der->test_function_type = order;
- pnode_der->name = varname;
- pnode_der->interpolate_name_der = pnode_der->interpolate_name;
- pnode_der->interpolate_name = interpolatename;
- }
+
+ case GA_NODE_X:
+ pnode->node_type = GA_NODE_CONSTANT;
+ pnode->init_matrix_tensor(meshdim, meshdim);
+ gmm::clear(pnode->tensor().as_vector());
+ for (size_type i = 0; i < meshdim; ++i)
+ pnode->tensor()(i,i) = scalar_type(1);
break;
+ case GA_NODE_NORMAL:
case GA_NODE_INTERPOLATE_NORMAL:
- GMM_ASSERT1(false, "Sorry, cannot derive the interpolated Normal");
- break;
+ GMM_ASSERT1(false, "Sorry, Gradient of Normal vector not implemented");
case GA_NODE_INTERPOLATE_DERIVATIVE:
- GMM_ASSERT1(false, "Sorry, second order transformation derivative "
- "not taken into account");
+ GMM_ASSERT1(false, "Sorry, gradient of the derivative of a "
+ "tranformation not implemented");
break;
case GA_NODE_INTERPOLATE_FILTER:
- ga_node_derivation(tree, workspace, m, child0, varname,
- interpolatename, order);
+ ga_node_grad(tree, workspace, m, child0);
break;
case GA_NODE_ELEMENTARY_VAL:
+ pnode->node_type = GA_NODE_ELEMENTARY_GRAD; break;
case GA_NODE_ELEMENTARY_GRAD:
+ pnode->node_type = GA_NODE_ELEMENTARY_HESS; break;
case GA_NODE_ELEMENTARY_HESS:
- case GA_NODE_ELEMENTARY_DIVERG:
+ GMM_ASSERT1(false, "Sorry, cannot derive an Hessian once more");
+ case GA_NODE_ELEMENTARY_DIVERG: // Hess_u : Id(meshdim)
+ pnode->node_type = GA_NODE_ELEMENTARY_HESS;
+ tree.duplicate_with_operation(pnode, GA_COLON);
+ child0 = pnode; pnode = pnode->parent; child1 = pnode->children[1];
+ child1->init_matrix_tensor(meshdim, meshdim);
+ gmm::clear(pnode->tensor().as_vector());
+ for (size_type i = 0; i < meshdim; ++i)
+ child1->tensor()(i,i) = scalar_type(1);
+ child1->node_type = GA_NODE_CONSTANT;
+ break;
+
case GA_NODE_XFEM_PLUS_VAL:
+ pnode->node_type = GA_NODE_XFEM_PLUS_GRAD; break;
case GA_NODE_XFEM_PLUS_GRAD:
+ pnode->node_type = GA_NODE_XFEM_PLUS_HESS; break;
case GA_NODE_XFEM_PLUS_HESS:
- case GA_NODE_XFEM_PLUS_DIVERG:
+ GMM_ASSERT1(false, "Sorry, cannot derive an Hessian once more");
+ case GA_NODE_XFEM_PLUS_DIVERG: // Hess_u : Id(meshdim)
+ pnode->node_type = GA_NODE_XFEM_PLUS_HESS;
+ tree.duplicate_with_operation(pnode, GA_COLON);
+ child0 = pnode; pnode = pnode->parent; child1 = pnode->children[1];
+ child1->init_matrix_tensor(meshdim, meshdim);
+ gmm::clear(pnode->tensor().as_vector());
+ for (size_type i = 0; i < meshdim; ++i)
+ child1->tensor()(i,i) = scalar_type(1);
+ child1->node_type = GA_NODE_CONSTANT;
+ break;
+
case GA_NODE_XFEM_MINUS_VAL:
+ pnode->node_type = GA_NODE_XFEM_MINUS_GRAD; break;
case GA_NODE_XFEM_MINUS_GRAD:
+ pnode->node_type = GA_NODE_XFEM_MINUS_HESS; break;
case GA_NODE_XFEM_MINUS_HESS:
- case GA_NODE_XFEM_MINUS_DIVERG:
- mi.resize(1); mi[0] = 2;
- for (size_type i = 0; i < pnode->tensor_order(); ++i)
- mi.push_back(pnode->tensor_proper_size(i));
- pnode->t.adjust_sizes(mi);
- switch(pnode->node_type) {
- case GA_NODE_ELEMENTARY_VAL:
- pnode->node_type = GA_NODE_ELEMENTARY_VAL_TEST; break;
- case GA_NODE_ELEMENTARY_GRAD:
- pnode->node_type = GA_NODE_ELEMENTARY_GRAD_TEST; break;
- case GA_NODE_ELEMENTARY_HESS:
- pnode->node_type = GA_NODE_ELEMENTARY_HESS_TEST; break;
- case GA_NODE_ELEMENTARY_DIVERG:
- pnode->node_type = GA_NODE_ELEMENTARY_DIVERG_TEST; break;
- case GA_NODE_XFEM_PLUS_VAL:
- pnode->node_type = GA_NODE_XFEM_PLUS_VAL_TEST; break;
- case GA_NODE_XFEM_PLUS_GRAD:
- pnode->node_type = GA_NODE_XFEM_PLUS_GRAD_TEST; break;
- case GA_NODE_XFEM_PLUS_HESS:
- pnode->node_type = GA_NODE_XFEM_PLUS_HESS_TEST; break;
- case GA_NODE_XFEM_PLUS_DIVERG:
- pnode->node_type = GA_NODE_XFEM_PLUS_DIVERG_TEST; break;
- case GA_NODE_XFEM_MINUS_VAL:
- pnode->node_type = GA_NODE_XFEM_MINUS_VAL_TEST; break;
- case GA_NODE_XFEM_MINUS_GRAD:
- pnode->node_type = GA_NODE_XFEM_MINUS_GRAD_TEST; break;
- case GA_NODE_XFEM_MINUS_HESS:
- pnode->node_type = GA_NODE_XFEM_MINUS_HESS_TEST; break;
- case GA_NODE_XFEM_MINUS_DIVERG:
- pnode->node_type = GA_NODE_XFEM_MINUS_DIVERG_TEST; break;
- default : GMM_ASSERT1(false, "internal error");
- }
- pnode->test_function_type = order;
+ GMM_ASSERT1(false, "Sorry, cannot derive an Hessian once more");
+ case GA_NODE_XFEM_MINUS_DIVERG: // Hess_u : Id(meshdim)
+ pnode->node_type = GA_NODE_XFEM_MINUS_HESS;
+ tree.duplicate_with_operation(pnode, GA_COLON);
+ child0 = pnode; pnode = pnode->parent; child1 = pnode->children[1];
+ child1->init_matrix_tensor(meshdim, meshdim);
+ gmm::clear(pnode->tensor().as_vector());
+ for (size_type i = 0; i < meshdim; ++i)
+ child1->tensor()(i,i) = scalar_type(1);
+ child1->node_type = GA_NODE_CONSTANT;
break;
case GA_NODE_OP:
switch(pnode->op_type) {
- case GA_PLUS: case GA_MINUS:case GA_NODE_DIVERG:
- if (mark0 && mark1) {
- ga_node_derivation(tree, workspace, m, child0, varname,
- interpolatename, order);
- ga_node_derivation(tree, workspace, m, child1, varname,
- interpolatename, order);
- } else if (mark0) {
- ga_node_derivation(tree, workspace, m, child0, varname,
- interpolatename, order);
- tree.replace_node_by_child(pnode, 0);
- } else {
- ga_node_derivation(tree, workspace, m, child1, varname,
- interpolatename, order);
- if (pnode->op_type == GA_MINUS) {
- pnode->op_type = GA_UNARY_MINUS;
- tree.clear_node(child0);
- }
- else
- tree.replace_node_by_child(pnode, 1);
- }
- break;
+ case GA_PLUS: case GA_MINUS:
+ if (mark0 && mark1) {
+ ga_node_grad(tree, workspace, m, child0);
+ ga_node_grad(tree, workspace, m, child1);
+ } else if (mark0) {
+ ga_node_grad(tree, workspace, m, child0);
+ tree.replace_node_by_child(pnode, 0);
+ } else {
+ ga_node_grad(tree, workspace, m, child1);
+ if (pnode->op_type == GA_MINUS) {
+ pnode->op_type = GA_UNARY_MINUS;
+ tree.clear_node(child0);
+ }
+ else
+ tree.replace_node_by_child(pnode, 1);
+ }
+ break;
+
+ case GA_UNARY_MINUS: case GA_PRINT:
+ ga_node_grad(tree, workspace, m, child0);
+ break;
- case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
- case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
- ga_node_derivation(tree, workspace, m, child0, varname,
+
+
+ #ifdef continue_here
+
+ case GA_QUOTE: case GA_SYM: case GA_SKEW:
+ case GA_TRACE: case GA_DEVIATOR:
+ ga_node_grad(tree, workspace, m, child0, varname,
interpolatename, order);
break;
@@ -3481,7 +3642,7 @@ namespace getfem {
if (mark0 && mark1) {
if (sub_tree_are_equal(child0, child1, workspace, 0) &&
(pnode->op_type != GA_MULT || child0->tensor_order() < 2)) {
- ga_node_derivation(tree, workspace, m, child1, varname,
+ ga_node_grad(tree, workspace, m, child1, varname,
interpolatename, order);
tree.insert_node(pnode, GA_NODE_OP);
pnode->parent->op_type = GA_MULT;
@@ -3496,17 +3657,17 @@ namespace getfem {
(child0->tensor_proper_size()== 1 &&
child1->tensor_proper_size()== 1))
std::swap(pnode->children[0], pnode->children[1]);
- ga_node_derivation(tree, workspace, m, child0, varname,
+ ga_node_grad(tree, workspace, m, child0, varname,
interpolatename, order);
- ga_node_derivation(tree, workspace, m,
+ ga_node_grad(tree, workspace, m,
pnode->parent->children[1]->children[1],
varname, interpolatename, order);
}
} else if (mark0) {
- ga_node_derivation(tree, workspace, m, child0, varname,
+ ga_node_grad(tree, workspace, m, child0, varname,
interpolatename, order);
} else
- ga_node_derivation(tree, workspace, m, child1, varname,
+ ga_node_grad(tree, workspace, m, child1, varname,
interpolatename, order);
break;
@@ -3518,7 +3679,7 @@ namespace getfem {
else {
if (mark0) {
tree.duplicate_with_subtraction(pnode);
- ga_node_derivation(tree, workspace, m, child0, varname,
+ ga_node_grad(tree, workspace, m, child0, varname,
interpolatename, order);
pnode = pnode->parent->children[1];
} else {
@@ -3541,22 +3702,23 @@ namespace getfem {
pnode_mult->children.resize(2, nullptr);
tree.copy_node(pnode_param->children[1],
pnode_mult, pnode_mult->children[1]);
- ga_node_derivation(tree, workspace, m, pnode_mult->children[1],
+ ga_node_grad(tree, workspace, m, pnode_mult->children[1],
varname, interpolatename, order);
} else {
- ga_node_derivation(tree, workspace, m, child0, varname,
+ ga_node_grad(tree, workspace, m, child0, varname,
interpolatename, order);
}
break;
-
+#endif
default: GMM_ASSERT1(false, "Unexpected operation. Internal error.");
}
break;
-
+
+#ifdef continue_here
case GA_NODE_C_MATRIX:
for (size_type i = 0; i < pnode->children.size(); ++i) {
if (pnode->children[i]->marked)
- ga_node_derivation(tree, workspace, m, pnode->children[i],
+ ga_node_grad(tree, workspace, m, pnode->children[i],
varname, interpolatename, order);
else {
pnode->children[i]->init_scalar_tensor(scalar_type(0));
@@ -3568,8 +3730,10 @@ namespace getfem {
case GA_NODE_PARAMS:
if (child0->node_type == GA_NODE_RESHAPE) {
- ga_node_derivation(tree, workspace, m, pnode->children[1],
+ ga_node_grad(tree, workspace, m, pnode->children[1],
varname, interpolatename, order);
+ } else if (child0->node_type == GA_NODE_CONTRACT) {
+ // TODO !!!! (avec mark1 et "child2"->marked
} else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
std::string name = child0->name;
ga_predef_function_tab::const_iterator it =
PREDEF_FUNCTIONS.find(name);
@@ -3631,7 +3795,7 @@ namespace getfem {
pnode_op->op_type = GA_DOTMULT;
pnode_op->children.resize(2, nullptr);
tree.copy_node(child1, pnode_op, pnode_op->children[1]);
- ga_node_derivation(tree, workspace, m, pnode_op->children[1],
+ ga_node_grad(tree, workspace, m, pnode_op->children[1],
varname, interpolatename, order);
}
} else {
@@ -3669,7 +3833,7 @@ namespace getfem {
pnode_op->op_type = GA_DOTMULT;
pnode_op->children.resize(2, nullptr);
tree.copy_node(child1, pnode_op, pnode_op->children[1]);
- ga_node_derivation(tree, workspace, m, pnode_op->children[1],
+ ga_node_grad(tree, workspace, m, pnode_op->children[1],
varname, interpolatename, order);
}
if (child2->marked) {
@@ -3706,7 +3870,7 @@ namespace getfem {
pnode_op->op_type = GA_DOTMULT;
pnode_op->children.resize(2, nullptr);
tree.copy_node(child2, pnode_op, pnode_op->children[1]);
- ga_node_derivation(tree, workspace, m, pnode_op->children[1],
+ ga_node_grad(tree, workspace, m, pnode_op->children[1],
varname, interpolatename, order);
}
}
@@ -3755,7 +3919,7 @@ namespace getfem {
pnode_op->children.resize(2, nullptr);
tree.copy_node(pnode->children[i], pnode_op,
pnode_op->children[1]);
- ga_node_derivation(tree, workspace, m, pnode_op->children[1],
+ ga_node_grad(tree, workspace, m, pnode_op->children[1],
varname, interpolatename, order);
if (pnode2->children[0]->name.compare("Norm_sqr") == 0
@@ -3890,6 +4054,19 @@ namespace getfem {
case GA_NODE_PARAMS:
if (child0->node_type == GA_NODE_RESHAPE)
return ga_node_is_affine(child1);
+ if (child0->node_type == GA_NODE_CONTRACT) {
+ if (pnode->children.size() == 4) {
+ return ga_node_is_affine(child1);
+ } else if (pnode->children.size() == 5) {
+ if (mark1 && pnode->children[3]->marked) return false;
+ if (mark1) return ga_node_is_affine(child1);
+ return ga_node_is_affine(pnode->children[3]);
+ } else if (pnode->children.size() == 7) {
+ if (mark1 && pnode->children[4]->marked) return false;
+ if (mark1) return ga_node_is_affine(child1);
+ return ga_node_is_affine(pnode->children[4]);
+ }
+ }
if (child0->node_type == GA_NODE_PREDEF_FUNC)
return false;
if (child0->node_type == GA_NODE_OPERATOR)
diff --git a/src/getfem_generic_assembly_tree.cc
b/src/getfem_generic_assembly_tree.cc
index d0bb3e5..157a6f8 100644
--- a/src/getfem_generic_assembly_tree.cc
+++ b/src/getfem_generic_assembly_tree.cc
@@ -1044,6 +1044,11 @@ namespace getfem {
GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
break;
+ case GA_NODE_CONTRACT:
+ str << "Contract";
+ GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
+ break;
+
case GA_NODE_C_MATRIX:
{
GMM_ASSERT1(pnode->children.size(), "Invalid tree");
@@ -1124,12 +1129,6 @@ namespace getfem {
// 3 : reserved prefix Dot and Previous
int ga_check_name_validity(const std::string &name) {
- if (!(name.compare("X")) ||
- !(name.compare("Normal")) ||
- !(name.compare("Reshape")))
- return 1;
-
-
if (name.compare(0, 11, "Derivative_") == 0)
return 2;