FullMatrix::copy_to only works for dim <= 3

[Tao Jin]

Dear all,

When I am implementing a spectrum decomposition for an arbitrary symmetric matrix, I need to copy a full matrix into a rank-2 tensor. I know that FullMatrix has a copy_to() member function like below

template<typename number >
template<int dim>
void FullMatrix< number >::copy_to ( Tensor< 2, dim > & T,
const size_type src_r_i = 0,
const size_type src_r_j = dim - 1,
const size_type src_c_i = 0,
const size_type src_c_j = dim - 1,
const unsigned int dst_r = 0,
const unsigned int dst_c = 0
) const

However, the above copy_to() function only works when the tensor dimension is less or equal to 3, that is,

Tensor<2, dim>, dim <=3.

When dim is larger than 3, I have to use for loops to copy each matrix entry into the rank-2 tensor.

I know that it is typical to set dim as 2 or 3. However, does it make sense to make FullMatrix::copy_to() work for a general rank-2 tensor?

By the way, in the file full_matrix.inst.in,

for (S : REAL_SCALARS)
  {
    template void FullMatrix<S>::copy_to<1>(Tensor<2, 1> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;

    template void FullMatrix<S>::copy_to<2>(Tensor<2, 2> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;

    template void FullMatrix<S>::copy_to<3>(Tensor<2, 3> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;
  }

This is probably the reason why copy_to only works for dim <=3

Best,

Tao

[Wolfgang Bangerth]

On 6/6/23 15:46, Tao Jin wrote:
>
> However, the above copy_to() function only works when the tensor dimension is
> less or equal to 3, that is,
> Tensor<2, dim>, dim <=3.
> When dim is larger than 3, I have to use for loops to copy each matrix entry
> into the rank-2 tensor.

What happens for dim>3? You only say "it only works when", but you don't
describe what goes wrong if dim>3?

> I know that it is typical to set dim as 2 or 3. However, does it make sense to
> make FullMatrix::copy_to() work for a general rank-2 tensor?

I see no reason why this should not work. The function is implemented in the
full_matrix.h file, so I believe it should just work if you call the function
on tensors that are larger than 3x3.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Tao Jin]

Dear Wolfgang,

Here is a simple test code (I also attached the source code):

int main()
{
  using namespace dealii;

  const unsigned int matrix_dimension = 4;

  Tensor<2, matrix_dimension> myTensor;

  FullMatrix<double> myMatrix(matrix_dimension,
      matrix_dimension);

  for (unsigned int i = 0; i < matrix_dimension; i++)
    for (unsigned int j = 0; j < matrix_dimension; j++)
      {
myMatrix(i,j) = i + j;
      }

  myMatrix.copy_to(myTensor);

  std::cout << myTensor.dimension << std::endl;

  return 0;
}

When matrix_dimension =2 or 3, the code works fine. When matrix_dimension = 4, it has a linking error as below:

error: undefined reference to 'void dealii::FullMatrix<double>::copy_to<4>(dealii::Tensor<2, 4, double>&, unsigned long, unsigned long, unsigned long, unsigned long, unsigned int, unsigned int) const'
collect2: error: ld returned 1 exit status
make[3]: *** [CMakeFiles/main.dir/build.make:116: main] Error 1

The version of deal.ii is 9.4.0.

I looked into the source code. Indeed, in "full_matrix.h", there is a declaration of "copy_to". However, the actual implementation is not in "/dealii-9.4.0/src/source/lac/full_matrix.cc", but in "/dealii-9.4.0/src/source/lac/full_matrix.inst.in"

for (S : REAL_SCALARS)
  {
    template void FullMatrix<S>::copy_to<1>(Tensor<2, 1> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;

    template void FullMatrix<S>::copy_to<2>(Tensor<2, 2> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;

    template void FullMatrix<S>::copy_to<3>(Tensor<2, 3> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;
  }

I assume this is related with explicit instantiation?

Best,

Tao

[Wolfgang Bangerth]

Tao:

On 6/7/23 12:45, Tao Jin wrote:
> *** Caution: EXTERNAL Sender ***

>
> Dear Wolfgang,
>
> Here is a simple test code (I also attached the source code):
>
> int main()
> {
>   using namespace dealii;
>

>   const unsigned int *matrix_dimension* = 4;

>
>   Tensor<2, matrix_dimension> myTensor;
>
>   FullMatrix<double> myMatrix(matrix_dimension,
>       matrix_dimension);
>
>   for (unsigned int i = 0; i < matrix_dimension; i++)
>     for (unsigned int j = 0; j < matrix_dimension; j++)
>       {
> myMatrix(i,j) = i + j;
>       }
>
>   myMatrix.copy_to(myTensor);
>
>   std::cout << myTensor.dimension << std::endl;
>
>   return 0;
> }
>

> When *matrix_dimension =*2 or 3, the code works fine. When *matrix_dimension
> *= 4, it has a linking error as below:
> /error: undefined reference to 'void

> dealii::FullMatrix<double>::copy_to<4>(dealii::Tensor<2, 4, double>&, unsigned
> long, unsigned long, unsigned long, unsigned long, unsigned int, unsigned int)
> const'
> collect2: error: ld returned 1 exit status

> make[3]: *** [CMakeFiles/main.dir/build.make:116: main] Error 1/

>
> The version of deal.ii is 9.4.0.

I tried this out with the current development version (which will become 9.5
within the next few weeks) and it works just fine:
https://github.com/dealii/dealii/pull/15331

I am not sure what is going on here. In your code snippet, you call the
copy_to() function with one argument, but the error message refers to the
function with 7 arguments. Are you sure that this is the code you wanted to use?

[Tao Jin]

Hi Wolfgang,

Below is the interface for copy_to:

template<typename number >
template<int dim>
void FullMatrix< number >::copy_to ( Tensor< 2, dim > & T,
const size_type src_r_i = 0,
const size_type src_r_j = dim - 1,
const size_type src_c_i = 0,
const size_type src_c_j = dim - 1,
const unsigned int dst_r = 0,
const unsigned int dst_c = 0
) const

When I called the function, I just used the default values for variables 2 to 7.

I actually compared the header file "include/deal.II/lac/full_matrix.h". In the develop branch (9.5.0), the definition of copy_to is implemented in the header file. However, in 9.4.0, the header file only has the declaration of copy_to. The actual implementation happens in "/dealii-9.4.0/src/source/lac/full_matrix.inst.in"

for (S : REAL_SCALARS)
  {
    template void FullMatrix<S>::copy_to<1>(Tensor<2, 1> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;

    template void FullMatrix<S>::copy_to<2>(Tensor<2, 2> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;

    template void FullMatrix<S>::copy_to<3>(Tensor<2, 3> &,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const size_type,
                                            const unsigned int,
                                            const unsigned int) const;
  }

This probably explains why the function works properly in 9.5.0 but does not work in 9.4.0 for dimensions larger than 3.

Anyhow, I think 9.5.0 has fixed this issue already.

Best,

Tao