Loading ...
Sorry, an error occurred while loading the content.

Re: [ublas-dev] ublas_compressed traits

Expand Messages
  • Gunter Winkler
    ... try adding #include or at-least the standard-definition of compressed matrix template
    Message 1 of 12 , Jul 21, 2003
    • 0 Attachment
      Henrik Dohlmann wrote:
      > #include <boost/numeric/bindings/traits/ublas_compressed.hpp>
      > namespace ublas = boost::numeric::ublas;
      >
      > int main(void)
      > {
      > ublas::compressed_matrix(6,6,18);
      > return 0;
      > }

      try adding
      #include <boost/numerik/ublas/config.hpp>

      or at-least the standard-definition of compressed matrix

      template<class T, class F = row_major, std::size_t IB = 0, class IA =
      unbounded_array<std::size_t>, class TA = unbounded_array<T> >
      class compressed_matrix;

      or use

      typedef compressed_matrix<double, row_major, 0,
      unbounded_array<std::size_t>, unbounded_array<double>
      atlas_compressed_matrix;

      HTH
      Gunter
    • Henrik Dohlmann
      Hi Gunter, ... Didn t help. There are still syntax errors in these lines from : //
      Message 2 of 12 , Jul 21, 2003
      • 0 Attachment
        Hi Gunter,

        >try adding
        >#include <boost/numeric/ublas/config.hpp>

        Didn't help. There are still syntax errors in these lines from
        <boost/numeric/bindings/traits/ublas_compressed.hpp>:

        // ublas::compressed_matrix<>
        template <typename T, typename F, std::size_t IB, typename IA, typename TA>
        *) struct sparse_matrix_traits<
        boost::numeric::ublas::compressed_matrix<T, F, IB, IA, TA>
        >

        At the point *) I get:
        boost\numeric\bindings\traits\ublas_compressed.hpp(32) : error C2143: syntax
        error : missing ';' before '<'
        boost\numeric\bindings\traits\ublas_compressed.hpp(32) : error C2059: syntax
        error : '<'


        The new and revised simple example is attached.

        Regards,
        Henrik
      • Henrik Dohlmann
        Hi Julius, ... Yeah, my example was a bit too quick :-) But the error is in the included header-file, so it didn t even reach it... Regards, Henrik
        Message 3 of 12 , Jul 21, 2003
        • 0 Attachment
          Hi Julius,

          >try
          >ublas::compressed_matrix<double>(6,6,18);
          >instead of
          >ublas::compressed_matrix(6,6,18);

          Yeah, my example was a bit too quick :-)
          But the error is in the included header-file, so it didn't even reach it...

          Regards,
          Henrik
        • kyriales
          Hi Henrik, you wrote: [...] ... Well, first of all, sparse_traits.hpp and ublas_compressed.hpp shuldn t be there (in bindings*.zip) at all. They are more than
          Message 4 of 12 , Jul 21, 2003
          • 0 Attachment
            Hi Henrik,

            you wrote:

            [...]
            > We had hoped that switching to compressed matrix instead would
            > reduce the time tremendously, but it seems we have stumbled upon
            > another strange Microsoft Compiler problem (MS Visual Studio .NET
            > 2003, aka VC++ 7.1).

            Well, first of all, sparse_traits.hpp and ublas_compressed.hpp
            shuldn't be there (in bindings*.zip) at all. They are more than
            experimental, part of my work on SuperLU bindings.

            Compressed matrices cannot be used with ATLAS (or BLAS
            or LAPACK), because ATLAS works only with dense matrices.

            SuperLU (http://crd.lbl.gov/~xiaoye/SuperLU/), on the other hand,
            works with compressed matrices.

            > Here is an extremely simple test-program to illustrate the trouble
            > we have:
            [...]
            > This results in the following errors:

            I think that the errors are caused by the fact that the inclusion of
            boost/numeric/bindings/traits/sparse_traits.hpp is missing in the
            ublas_compressed.hpp. (Thanx for cathing the bug ;o).

            Regards,

            fres
          • kyriales
            Hi again, (unrelated to ublas_compressed traits problem) ... bindings, but ... with a ... Why do you invert stiffness matrix? Numerical analysis books
            Message 5 of 12 , Jul 21, 2003
            • 0 Attachment
              Hi again,

              (unrelated to ublas_compressed traits problem)

              > We are working hard on a FEM implementation using uBlas and Atlas
              bindings, but
              > currently our static solution takes around 30 seconds for a problem
              with a
              > matrix of approximately dimension 2300x2300.
              >
              > ublas::vector<int> ipiv (nvars); // pivot vector
              > int rc = atlas::lu_factor (Ktmp, ipiv); // alias for getrf()
              > assert(!rc);
              > atlas::lu_invert(Ktmp, ipiv);
              > atlas::gemv(Ktmp,f,u);

              Why do you invert stiffness matrix? Numerical analysis books
              recommend back substitution:

              int rc = atlas::lu_factor (Ktmp, ipiv);
              assert (rc);
              // copy vector f to first (and only) column of matrix u
              atlas::lu_substitute (Ktmp, ipiv, u); // alias for getrs()

              It is more numerically stable and probably faster.

              Regards,

              fres
            • Henrik Dohlmann
              Hi, ... Okay, we tried that. But our application dies horribly when it enters atlas::lu_substitute. Any idea why? According to doc, our U has to be (nrhs by
              Message 6 of 12 , Jul 22, 2003
              • 0 Attachment
                Hi,

                >Why do you invert stiffness matrix? Numerical analysis books
                >recommend back substitution:
                >
                > int rc = atlas::lu_factor (Ktmp, ipiv);
                > assert (rc);
                > // copy vector f to first (and only) column of matrix u
                > atlas::lu_substitute (Ktmp, ipiv, u); // alias for getrs()
                >
                >It is more numerically stable and probably faster.

                Okay, we tried that. But our application dies horribly when it enters
                atlas::lu_substitute.
                Any idea why? According to doc, our U has to be (nrhs by n), right?

                int rc = atlas::lu_factor (Ktmp, ipiv); // alias for getrf()
                assert (!rc);
                matrix<scalar, row_major> U (1,nvars);
                row (U,1) = f;
                rc = atlas::lu_substitute (Ktmp, ipiv, U); // alias for getrs()
                u = row(U,1);

                Regards,
                Henrik
              • Kresimir Fresl
                Hi Henrik, ... [...] ... Yes. ... row (U, 0) = f; ... u = row (U, 0); But garbage in U can t be a reason for getrs() to die. How does it die? Error message?
                Message 7 of 12 , Jul 22, 2003
                • 0 Attachment
                  Hi Henrik,

                  >>Why do you invert stiffness matrix? Numerical analysis books
                  >>recommend back substitution:
                  [...]
                  > Okay, we tried that. But our application dies horribly when it enters
                  > atlas::lu_substitute.
                  > Any idea why? According to doc, our U has to be (nrhs by n), right?

                  Yes.

                  > int rc = atlas::lu_factor (Ktmp, ipiv); // alias for getrf()
                  > assert (!rc);
                  > matrix<scalar, row_major> U (1,nvars);
                  > row (U,1) = f;

                  row (U, 0) = f;

                  > rc = atlas::lu_substitute (Ktmp, ipiv, U); // alias for getrs()
                  > u = row(U,1);

                  u = row (U, 0);

                  But garbage in U can't be a reason for getrs() to die.
                  How does it die? Error message? Program hangs?
                  Program aborts? Or you simply get garbage?

                  Regards,

                  fres
                • Henrik Dohlmann
                  Hi, ... Aargh, I think the brand of coffee we are drinking, doen t work as expected :-) ... It just aborts! For some weird reason, we cannot get VC71 to enter
                  Message 8 of 12 , Jul 22, 2003
                  • 0 Attachment
                    Hi,

                    >row (U, 0) = f;

                    Aargh, I think the brand of coffee we are drinking, doen't work as expected :-)

                    >But garbage in U can't be a reason for getrs() to die.
                    >How does it die? Error message? Program hangs?
                    >Program aborts? Or you simply get garbage?

                    It just aborts! For some weird reason, we cannot get VC71 to enter debug mode
                    anymore. It says that the generated symbol file is corrupt. And we have no clue
                    why?

                    >Why do you invert stiffness matrix?
                    >Numerical analysis books recommend back substitution:
                    >It is more numerically stable and probably faster.

                    It gave a tremendous speed-up. From 18 seconds downto 5 seconds.

                    Thanks.

                    Now, we "just" need to drop atlas and go for a sparse solution instead.

                    What is the big difference between compressed_matrix and sparse_matrix anyhow?

                    Regards,
                    Henrik
                  • Kresimir Fresl
                    Hi again, ... Maybe you should try to put some cognac or whisky in your coffee ;o) ... When you had `row (U, 1) = f; part of the memory which maybe contained
                    Message 9 of 12 , Jul 22, 2003
                    • 0 Attachment
                      Hi again,

                      >>row (U, 0) = f;

                      > Aargh, I think the brand of coffee we are drinking, doen't work as expected :-)

                      Maybe you should try to put some cognac or whisky in your coffee ;o)

                      >>But garbage in U can't be a reason for getrs() to die.
                      >>How does it die? Error message? Program hangs?
                      >>Program aborts? Or you simply get garbage?

                      > It just aborts! For some weird reason, we cannot get VC71 to enter debug mode
                      > anymore. It says that the generated symbol file is corrupt. And we have no clue
                      > why?

                      When you had `row (U, 1) = f;' part of the memory which maybe
                      contained program code was overwritten.

                      >>Why do you invert stiffness matrix?
                      >>Numerical analysis books recommend back substitution:
                      >>It is more numerically stable and probably faster.

                      > It gave a tremendous speed-up. From 18 seconds downto 5 seconds.

                      Your stiffness matrix is probably symmetric positive definite.
                      So you can try potrf()/potrs().

                      > Thanks.

                      You're welcome ;o)

                      > Now, we "just" need to drop atlas and go for a sparse solution instead.
                      >
                      > What is the big difference between compressed_matrix and sparse_matrix anyhow?

                      sparse_matrix uses some kind of associative array to store
                      elements (e.g. ublas::map_array or std::map).
                      compressed_matrix uses three arrays -- for row major storage
                      `pointers' (sorted) to first non-zero element in each row,
                      column indices (sorted) of non-zero elements in each row,
                      non-zero elements.

                      sparse_matrix with std::map allows O(log nnz) random insertations
                      (where nnz is number of non-zero elements). Insertations in other
                      matrix types require moving all elements which are after the insert
                      position. So, for e.g. global stiffness matrix assembly sparse_matrix
                      with std::map is `ideal' data structure -- assembly is O(nnz * log nnz)
                      operation. (There is also coordinate_matrix. It also uses three arrays:
                      row and column indices of non-zero elements and non-zero elements.
                      When elements are inserted, they are in fact appended at the end (O(1)
                      operation), so matrix become unsorted. But before other uses, matrix
                      must be sorted -- again O(nnz * log nnz) operation.)

                      On the other hand, matrix traversals are usually faster for
                      compressed matrix -- storage is contiguous and there's no
                      need to `chase' pointer. And, probably even more important,
                      Fortran and C libraries for solving linear systems with sparse
                      matrices (e.g. SuperLU) usually expect compressed (or, sometimes,
                      coordinate) storage formats.

                      So, you can use sparse_matrix<std::map> or coordinate_matrix
                      for stiffness matrix assembly and then copy it to compressed_matrix
                      to solve the system.

                      You can find more information about sparse matrices in excellent
                      book by Yousef Saad `Iterative methods for sparse linear systems'
                      -- free .ps and .pdf can be found at

                      http://www-users.cs.umn.edu/~saad/

                      (Unfortunately, as Saad uses Fortran, associative array approach
                      is not mentioned -- not that associative arrays are impossible in
                      Fortran, but it's not simple to implement them; and there are
                      other, equally efficient and more `natural' solutions.)

                      Regards,

                      fres
                    • Henrik Dohlmann
                      Hi, ... Right now it is, yes. This reduces computation time to just over 3 seconds: // This solves a symmetric version using Cholesky-factorization and
                      Message 10 of 12 , Jul 24, 2003
                      • 0 Attachment
                        Hi,

                        >> It gave a tremendous speed-up. From 18 seconds downto 5 seconds.
                        >Your stiffness matrix is probably symmetric positive definite.
                        >So you can try potrf()/potrs().

                        Right now it is, yes. This reduces computation time to just over 3 seconds:

                        // This solves a symmetric version using Cholesky-factorization and
                        back-substitution.
                        //--- See: ...\boost\libs\numeric\bindings\atlas\doc\index.html
                        matrix<scalar, row_major> U (1,nvars); // Matrix of size (nrhs x n)
                        symmetric_adaptor<matrix<scalar>, lower> SA(Ktmp); // Symmetric version of K
                        atlas::cholesky_factor (SA); // alias for potrf()
                        row (U,0) = f; // fill in "all" the rhs's
                        atlas::cholesky_substitute (SA, U); // alias for potrs()
                        u = row(U,0); // grab the solution

                        Thanks.

                        Regards,
                        Henrik
                      Your message has been successfully submitted and would be delivered to recipients shortly.