Hi, These can very well be very naive questions, since I am relatively new to this - apologies..: I am trying to write a pde solver for a 3D problem. Multi-array seems like a very natural choice for a data structure to store the functtion values. The method I am trying to use is operrator splitting. This amounts to solving a tridiagonal system for each dimension in sequence. Suppose then that I have a MA( M, N, L ): It is highly desirable to maintain this xyz (let's call it) "logical" system of reference since incorporation of boundary conditions and stencil calculations would be conceptually straightforward -otherwise very complicated. The idea being that, I would then have the tridiagonal opertator ( and rhs) in the "rotated" system and all other calculations in the "logical" one. Therefore, I need some fascility to rotate my MA, from the xyz to yzx etc (cyclically). Is there functionality ( as a view to the data ) that would allow me to do this off-the-shelve? Also, a couple of other points: -If I wanted to use views to subgrids, in parallelized calculations, are there memory allocations that happen under the hood for which I should worry (for performance)? -Is a one dimensional MA usable? -Are there adapters between 2D MAs and uBlas matrices? Thank you very much for your help in advance! P-
Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form. I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros -----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes On 05/21/11 15:02, petros wrote: Hi Petros, [snip]
Suppose then that I have a MA( M, N, L ):
Please be more specific. I assume: MA means boost::multi_array however, I've no idea what M, N, and L are. Could you please explain further? -regards, Larry _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hello Petros, I think that what you are looking for is a different storage order. Take a look at "boost/multi_array/storage_order.hpp". 3 classes are defined: general_storage_order, c_storage_order and fortran_storage_order. You need to make a "cyclic_storage_order" (or something like that) that receives an integer. When receiving zero, it should behave like c_storage_order. When receiving 1, everything should be permuted by 1 etc. This new class should inherit from general_storage_order. Use c_storage_order and fortran_storage_order as examples. Keep "ascending" true, only "ordering" has to be changed. A modulo should do it. Hope it helps, Pierre-André Noël On Sun, May 22, 2011 at 11:53 AM, petros <pmamales@nyc.rr.com> wrote:
Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form.
I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros
-----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes
On 05/21/11 15:02, petros wrote: Hi Petros,
[snip]
Suppose then that I have a MA( M, N, L ):
Please be more specific. I assume:
MA means boost::multi_array
however, I've no idea what M, N, and L are. Could you please explain further?
-regards, Larry
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
I forgot a detail: once you got your "cyclic_storage_order", you should make a normal multi_array as usual then access it through a multi_array_ref that uses "cyclic_storage_order". On Sun, May 22, 2011 at 12:17 PM, Pierre-Andre Noel < noel.pierre.andre@gmail.com> wrote:
Hello Petros,
I think that what you are looking for is a different storage order.
Take a look at "boost/multi_array/storage_order.hpp". 3 classes are defined: general_storage_order, c_storage_order and fortran_storage_order.
You need to make a "cyclic_storage_order" (or something like that) that receives an integer. When receiving zero, it should behave like c_storage_order. When receiving 1, everything should be permuted by 1 etc. This new class should inherit from general_storage_order. Use c_storage_order and fortran_storage_order as examples.
Keep "ascending" true, only "ordering" has to be changed. A modulo should do it.
Hope it helps,
Pierre-André Noël
On Sun, May 22, 2011 at 11:53 AM, petros <pmamales@nyc.rr.com> wrote:
Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form.
I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros
-----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes
On 05/21/11 15:02, petros wrote: Hi Petros,
[snip]
Suppose then that I have a MA( M, N, L ):
Please be more specific. I assume:
MA means boost::multi_array
however, I've no idea what M, N, and L are. Could you please explain further?
-regards, Larry
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hi Pierre-Andre, Thank you very much. Yes this helps immensely for the right hand-side of the derived system (say Ax = b). But for the creation of the matrix A (which incorporates boundary conditions and stencil calcs) a view with similar properties would still be needed. Also, can the starting point be a fortran storage order? Thank you again, Petros From: Pierre-Andre Noel Sent: Sunday, May 22, 2011 12:27 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes I forgot a detail: once you got your "cyclic_storage_order", you should make a normal multi_array as usual then access it through a multi_array_ref that uses "cyclic_storage_order". On Sun, May 22, 2011 at 12:17 PM, Pierre-Andre Noel <noel.pierre.andre@gmail.com> wrote: Hello Petros, I think that what you are looking for is a different storage order. Take a look at "boost/multi_array/storage_order.hpp". 3 classes are defined: general_storage_order, c_storage_order and fortran_storage_order. You need to make a "cyclic_storage_order" (or something like that) that receives an integer. When receiving zero, it should behave like c_storage_order. When receiving 1, everything should be permuted by 1 etc. This new class should inherit from general_storage_order. Use c_storage_order and fortran_storage_order as examples. Keep "ascending" true, only "ordering" has to be changed. A modulo should do it. Hope it helps, Pierre-André Noël On Sun, May 22, 2011 at 11:53 AM, petros <pmamales@nyc.rr.com> wrote: Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form. I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros -----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes On 05/21/11 15:02, petros wrote: Hi Petros, [snip] Suppose then that I have a MA( M, N, L ): Please be more specific. I assume: MA means boost::multi_array however, I've no idea what M, N, and L are. Could you please explain further? -regards, Larry _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users -------------------------------------------------------------------------------- _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hi Petros, Yes, you could use fortran_storage_order as a starting point. I'm not sure I understand the problem. Can't you use array_ref as your views? Here is an example (not tested), given the hypothetical "cyclic_storage_order". boost::multi_array<value_type,3> A(boost::extents[L][M][N]); boost::multi_array_ref<value_type,3> A_c1(A.data(), boost::extents[M][N][L],cyclic_storage_order(1)); boost::multi_array_ref<value_type,3> A_c2(A.data(), boost::extents[N][L][M],cyclic_storage_order(2)); Now all A, A_c1 and A_c2 use the same memory (and you can thus initialize A through A_c1 or A_c2). If you need a "rotated 2D view", do as usual but instead use A_c1 or A_c2. I hope this answers your question. Sincerely, Pierre-André Noël On Sun, May 22, 2011 at 1:34 PM, petros <pmamales@nyc.rr.com> wrote:
Hi Pierre-Andre, Thank you very much. Yes this helps immensely for the right hand-side of the derived system (say Ax = b). But for the creation of the matrix A (which incorporates boundary conditions and stencil calcs) a view with similar properties would still be needed. Also, can the starting point be a fortran storage order? Thank you again, Petros
*From:* Pierre-Andre Noel <noel.pierre.andre@gmail.com> *Sent:* Sunday, May 22, 2011 12:27 PM *To:* boost-users@lists.boost.org *Subject:* Re: [Boost-users] multi-array and pdes
I forgot a detail: once you got your "cyclic_storage_order", you should make a normal multi_array as usual then access it through a multi_array_ref that uses "cyclic_storage_order".
On Sun, May 22, 2011 at 12:17 PM, Pierre-Andre Noel < noel.pierre.andre@gmail.com> wrote:
Hello Petros,
I think that what you are looking for is a different storage order.
Take a look at "boost/multi_array/storage_order.hpp". 3 classes are defined: general_storage_order, c_storage_order and fortran_storage_order.
You need to make a "cyclic_storage_order" (or something like that) that receives an integer. When receiving zero, it should behave like c_storage_order. When receiving 1, everything should be permuted by 1 etc. This new class should inherit from general_storage_order. Use c_storage_order and fortran_storage_order as examples.
Keep "ascending" true, only "ordering" has to be changed. A modulo should do it.
Hope it helps,
Pierre-André Noël
On Sun, May 22, 2011 at 11:53 AM, petros <pmamales@nyc.rr.com> wrote:
Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form.
I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros
-----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes
On 05/21/11 15:02, petros wrote: Hi Petros,
[snip]
Suppose then that I have a MA( M, N, L ):
Please be more specific. I assume:
MA means boost::multi_array
however, I've no idea what M, N, and L are. Could you please explain further?
-regards, Larry
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
------------------------------ _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hi Pierre-Andre, This looks very promising and could save me a lot of work..). Will look into it. Thank you very much, Petros From: Pierre-Andre Noel Sent: Sunday, May 22, 2011 2:06 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes Hi Petros, Yes, you could use fortran_storage_order as a starting point. I'm not sure I understand the problem. Can't you use array_ref as your views? Here is an example (not tested), given the hypothetical "cyclic_storage_order". boost::multi_array<value_type,3> A(boost::extents[L][M][N]); boost::multi_array_ref<value_type,3> A_c1(A.data(), boost::extents[M][N][L],cyclic_storage_order(1)); boost::multi_array_ref<value_type,3> A_c2(A.data(), boost::extents[N][L][M],cyclic_storage_order(2)); Now all A, A_c1 and A_c2 use the same memory (and you can thus initialize A through A_c1 or A_c2). If you need a "rotated 2D view", do as usual but instead use A_c1 or A_c2. I hope this answers your question. Sincerely, Pierre-André Noël On Sun, May 22, 2011 at 1:34 PM, petros <pmamales@nyc.rr.com> wrote: Hi Pierre-Andre, Thank you very much. Yes this helps immensely for the right hand-side of the derived system (say Ax = b). But for the creation of the matrix A (which incorporates boundary conditions and stencil calcs) a view with similar properties would still be needed. Also, can the starting point be a fortran storage order? Thank you again, Petros From: Pierre-Andre Noel Sent: Sunday, May 22, 2011 12:27 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes I forgot a detail: once you got your "cyclic_storage_order", you should make a normal multi_array as usual then access it through a multi_array_ref that uses "cyclic_storage_order". On Sun, May 22, 2011 at 12:17 PM, Pierre-Andre Noel <noel.pierre.andre@gmail.com> wrote: Hello Petros, I think that what you are looking for is a different storage order. Take a look at "boost/multi_array/storage_order.hpp". 3 classes are defined: general_storage_order, c_storage_order and fortran_storage_order. You need to make a "cyclic_storage_order" (or something like that) that receives an integer. When receiving zero, it should behave like c_storage_order. When receiving 1, everything should be permuted by 1 etc. This new class should inherit from general_storage_order. Use c_storage_order and fortran_storage_order as examples. Keep "ascending" true, only "ordering" has to be changed. A modulo should do it. Hope it helps, Pierre-André Noël On Sun, May 22, 2011 at 11:53 AM, petros <pmamales@nyc.rr.com> wrote: Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form. I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros -----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes On 05/21/11 15:02, petros wrote: Hi Petros, [snip] Suppose then that I have a MA( M, N, L ): Please be more specific. I assume: MA means boost::multi_array however, I've no idea what M, N, and L are. Could you please explain further? -regards, Larry _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users ------------------------------------------------------------------------------ _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users -------------------------------------------------------------------------------- _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hi Pierre-Andre, Thank you for your response. For a moment, I thought that the cyclic_storage_order was already avail.. A question: why does it have to derive for the general_storage class? there is the implicit conversion operator, no? I did something like: class cyclic_storage_order { typedef detail::multi_array::size_type size_type; public: cyclic_storage_order( const std::size_t shift ) :shift_(shift) {} template <std::size_t NumDims> operator general_storage_order<NumDims>() const { boost::array<size_type,NumDims> ordering; boost::array<bool,NumDims> ascending; for (size_type i=0; i != NumDims; ++i) { ordering[i] = (i+shift_)%NumDims; ascending[i] = true; } return general_storage_order<NumDims>(ordering.begin(), ascending.begin()); } private: std::size_t shift_; }; Is this what you meant? Thank you very much for all your help, P- ---- Pierre-Andre Noel <noel.pierre.andre@gmail.com> wrote:
I forgot a detail: once you got your "cyclic_storage_order", you should make a normal multi_array as usual then access it through a multi_array_ref that uses "cyclic_storage_order".
On Sun, May 22, 2011 at 12:17 PM, Pierre-Andre Noel < noel.pierre.andre@gmail.com> wrote:
Hello Petros,
I think that what you are looking for is a different storage order.
Take a look at "boost/multi_array/storage_order.hpp". 3 classes are defined: general_storage_order, c_storage_order and fortran_storage_order.
You need to make a "cyclic_storage_order" (or something like that) that receives an integer. When receiving zero, it should behave like c_storage_order. When receiving 1, everything should be permuted by 1 etc. This new class should inherit from general_storage_order. Use c_storage_order and fortran_storage_order as examples.
Keep "ascending" true, only "ordering" has to be changed. A modulo should do it.
Hope it helps,
Pierre-André Noël
On Sun, May 22, 2011 at 11:53 AM, petros <pmamales@nyc.rr.com> wrote:
Hi Larry, apologies, for not being clear. Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet). I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists. My problem amounts to having a view that would allow for the permutation (if not rotation) of axes. As I begun to explain, the operator splitting methods decompose a problem in N dimensions into N 1-dimensional problems. The problem itself is stated, though, in a certain reference system (this means that the discretized operators (stencils) and the boundary conditions are expressed in terms of these boundary conditions. However, the linear equation systems that are produced, are in tridiagonal form, only in one of these steps ( the first step of the 1-D problems ). Otherwise it produces a block-banded system of equations, that can be reduced back to the familiar tridiagonal ( conceptually), through axes permutation (or rotation). If one things in terms of the 2D problem analog, while conceptually it would amount to exchange the two axes, in matrix form, and because of the storage alignment it would involve a much bigger operation. Since, these operators/systems are created from the problem, it would be very nice to have this view functionality to exploit for keeping the storage of the operators in the desired form.
I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, Petros
-----Original Message----- From: Larry Evans Sent: Sunday, May 22, 2011 10:43 AM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes
On 05/21/11 15:02, petros wrote: Hi Petros,
[snip]
Suppose then that I have a MA( M, N, L ):
Please be more specific. I assume:
MA means boost::multi_array
however, I've no idea what M, N, and L are. Could you please explain further?
-regards, Larry
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
On 05/22/11 10:53, petros wrote:
Hi Larry, apologies, for not being clear. No problem.
Yes, MA is boost::multi_array. M, N, L are the dimensions (this is only a "logical" statement, rather than a proper code snippet).
I think I have pretty much sorted out the calling practices wrt multi-aray and looking for if certain functionality exists.
My problem amounts to having a view that would allow for the permutation (if not rotation) of axes.
I'm not sure I understand permutation and rotation. For example, say you have: multi_array<T,2> ma2(extents[3][4]); The array would have 3 rows and 4 columns. Now, since there are only 2 axes (row and column) permuting the axes would be transposing the matrix. Is that right. In general, for: multi_array<T,N> maN(extents[L[0]][L[1]]...[L[(N-1)]]); a permutaion of the axes could be specified by a permutation, p, of: (0,1,...,N-1) and the resulting array, maT, would satisify, for all: unsigned i[N]; such that for all axis in (0...N-1), i[axis]<L[axis]: maN[i[0]][i[1]]...[i[N-1]] == maT[i[p[0]]][i[p[1]]]...[i[p[N-1]]] If this is what you mean, then this could be done using the stepper data structure mentioned here: http://article.gmane.org/gmane.comp.lib.boost.devel/218623/match=index_list The q_i in that post corresponds to the p vector (which is a permutation of (0...N-1) mentioned above. The book referenced in that post also has a section on rotation; hence, it may be useful for that also. [snip]
I hope this is clearer now - don't mean to drag you into my problem. Thank you for your response, [snip]
You're welcome, Petros. HTH. -Larry
On 05/22/11 12:17, Larry Evans wrote: [snip]
The book referenced in that post also has a section on rotation; hence, it may be useful for that also.
Just to be clear about rotation, suppose: multi_array<int,2> ma2(extents[3][4]); and suppose it was filled with 1...12 so that a printout would show: { { 1, 2, 3, 4} , { 5, 6, 7, 8} , { 9,10,11,12} } Then, to rotate the 3 rows by 2 1 and -1 elements would produce: { { 3, 4, 1, 2} , { 6, 7, 8, 5} , {12, 9,10,11} } and to rotate the 4 columns by 0 1 -1 -2 would produce: { { 1, 6,11, 8} , { 5,10, 3,12} , { 9, 2, 7, 4} } Is that what you meant by rotation? -regards, Larry
Hi Larry, Thank you very much for your extended response. I am not sure, will have to think about it, altough this seems right. In appreciation of your effort to help me, let me give you some color: Say I am trying to sove a 3d problem using splitting methods. Lets say that the original system f reference is xyz. One alays ends up to a system of equations in the vectorized reprezentation of the grid (very much like the array where the elements of the ma are stored). Then, when trying to solve the problem in the x direction (while in fortran storage scheme), I obtain a nice tridiagonal system of equations which I can solve very efficiently (using Thomas algorithm which is O(N) ). When I go to the second dimension, the tridiagonal system is hidden (in the original vector). However, in the rotsted yzx system it is there!! So, I am "paying the price" of writting it in the rotated view and copying it into a different vector to solve efficiently (o/wise I end up with a much more expensive sparse system to solve..) Hence the root of my questions. One of the issues is that the rows that I need to permute is rather complicated and with a little of smoke and mirrors, I hope to bypass the complications. Apologies, if all this is known to ou already or if went way too atray. Also, apologies for not catching your response earlier. Thank you agian for all your help-will look inth the references, Petros ---- Larry Evans <cppljevans@suddenlink.net> wrote:
On 05/22/11 12:17, Larry Evans wrote: [snip]
The book referenced in that post also has a section on rotation; hence, it may be useful for that also.
Just to be clear about rotation, suppose:
multi_array<int,2> ma2(extents[3][4]);
and suppose it was filled with 1...12 so that a printout would show:
{ { 1, 2, 3, 4} , { 5, 6, 7, 8} , { 9,10,11,12} }
Then, to rotate the 3 rows by 2 1 and -1 elements would produce:
{ { 3, 4, 1, 2} , { 6, 7, 8, 5} , {12, 9,10,11} }
and to rotate the 4 columns by 0 1 -1 -2 would produce:
{ { 1, 6,11, 8} , { 5,10, 3,12} , { 9, 2, 7, 4} }
Is that what you meant by rotation?
-regards, Larry
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Hi Larry,
Thank you very much for your extended response. You are very welcome! [snip] I went ahead an implemented a successor to the array_dyn mentioned before. The revised array_dyn is now in the boost vault under
On 05/24/11 00:58, pmamales@nyc.rr.com wrote: Hi Petros, the data structures directory: http://www.boostpro.com/vault/index.php?&directory=Data%20Structures It has one feature which you're probably not interested in, the axis reversal. To reverse an axis, you specify the length as a negative number. To rotate all the axes( the feature you are interested in), you use the axis_rot(rotation) method. The storage order is specified in with a CTOR arg of type enum dirs {dir_fwd, dir_rev}. I think these can be used to chose either fortran or c-like storage order. Hopefully the printout will help you decide which value of the parameter you desire. The test driver, in array_dyn.cpp, exercises these features and shows the results with printouts. The printout is shown in array_dyn.out member of the zip file. HTH. -Larry
Hi Larry,
Thank you very much for your extended response. You are very welcome! [snip] I went ahead an implemented a successor to the array_dyn mentioned before. The revised array_dyn is now in the boost vault under
Hi Larry, This is a lot of help! Thank you very very much! Will look into it. Best Regards, Petros -----Original Message----- From: Larry Evans Sent: Tuesday, May 24, 2011 3:54 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes On 05/24/11 00:58, pmamales@nyc.rr.com wrote: Hi Petros, the data structures directory: http://www.boostpro.com/vault/index.php?&directory=Data%20Structures It has one feature which you're probably not interested in, the axis reversal. To reverse an axis, you specify the length as a negative number. To rotate all the axes( the feature you are interested in), you use the axis_rot(rotation) method. The storage order is specified in with a CTOR arg of type enum dirs {dir_fwd, dir_rev}. I think these can be used to chose either fortran or c-like storage order. Hopefully the printout will help you decide which value of the parameter you desire. The test driver, in array_dyn.cpp, exercises these features and shows the results with printouts. The printout is shown in array_dyn.out member of the zip file. HTH. -Larry _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
On 05/24/11 16:09, petros wrote:
Hi Larry,
Hi Larry,
Thank you very much for your extended response. You are very welcome! [snip] I went ahead an implemented a successor to the array_dyn mentioned before. The revised array_dyn is now in the boost vault under
This is a lot of help! Thank you very very much! Will look into it. Best Regards, Petros -----Original Message----- From: Larry Evans Sent: Tuesday, May 24, 2011 3:54 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes On 05/24/11 00:58, pmamales@nyc.rr.com wrote: Hi Petros, the data structures directory:
http://www.boostpro.com/vault/index.php?&directory=Data%20Structures
It has one feature which you're probably not interested in, the
The attached test driver, when used with the vault code mentioned above, shows how the multi_array storage order compares with the enum dirs{dir_fwd,dir_rev} arg to array_dyn. The comments at the top summarize the comparison. The output showing the printouts is also attached. Maybe one advantage of array_dyn is the ease of rotating the axes. Another may be offset calculation, because with array_dyn, no testing is done for the bool values of general_storage_order::ascending. Instead, an offset and negative strides achieve the same result, AFAICT. -regards, Larry
On 05/25/11 15:44, Larry Evans wrote: [snip]
Maybe one advantage of array_dyn is the ease of rotating the axes. Another may be offset calculation, because with array_dyn, no testing is done for the bool values of general_storage_order::ascending. Instead, an offset and negative strides achieve the same result, AFAICT.
After finding: Reference access(boost::type<Reference>,index idx,TPtr base, const size_type* extents, const index* strides, const index* index_bases) const { BOOST_ASSERT(idx - index_bases[0] >= 0); BOOST_ASSERT(size_type(idx - index_bases[0]) < extents[0]); // return a sub_array<T,NDims-1> proxy object TPtr newbase = base + idx * strides[0]; return Reference(newbase,extents+1,strides+1,index_bases+1); in boost/multi_array/base.hpp, I'm guessing that array_dyn has no advantage in offset calculation since there's no test in that base.hpp code on general_storage_order<>::ascending. However, I'm still not sure how multi_array does it, unless it allows strides to be negative to account for a false ascending value and increments the base value to account for the negative strides, much like array_dyn's: box_domain::init_iter_strides Sorry for the noise about offset calculations. OTOH, there would be no need, in array_dyn, for incrementing all of extents, strides, and index_bases pointers since they're all (well, just extents and strides or really lengths and strides) stored in the same element of an array: box_domain<>::my_length_strides; ; hence, only one increment would be needed. Of course that's a small advantage, and there's probably disadvantages to this array_dyn which I've not thought of yet.
On 05/27/11 14:14, Larry Evans wrote: [snip]
After finding:
Reference access(boost::type<Reference>,index idx,TPtr base, const size_type* extents, const index* strides, const index* index_bases) const {
BOOST_ASSERT(idx - index_bases[0] >= 0); BOOST_ASSERT(size_type(idx - index_bases[0]) < extents[0]); // return a sub_array<T,NDims-1> proxy object TPtr newbase = base + idx * strides[0]; return Reference(newbase,extents+1,strides+1,index_bases+1);
in boost/multi_array/base.hpp, I'm guessing that array_dyn has no advantage in offset calculation since there's no test in that base.hpp code on general_storage_order<>::ascending. However, I'm still not sure how multi_array does it, unless it allows strides to be negative to account for a false ascending value and increments the base value to account for the negative strides, much like array_dyn's:
box_domain::init_iter_strides
Sorry for the noise about offset calculations.
After reading: http://www.boost.org/doc/libs/1_46_1/libs/multi_array /doc/reference.html#memory_layout which provides an example: Finally, both dimensions could be stored in descending order: int data[] = {11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0}; int *a = data + 11; int s[] = { -4, -1 }; which strongly suggest that multi_array *does* use negative strides: int s[] = { -4, -1 }; *and* somehow performs an offset calculation to account for those: int *a = data + 11; I'm guessing that there's really no difference in the way array_dyn implements this reversal and offset calculation. The only difference is in the interface which uses a bool to specify whether an axis is reversed rather than a negative length. I just didn't see where this was done from looking at code :( -regards Larry
On 05/24/11 00:58, pmamales@nyc.rr.com wrote:
Hi Larry,
Thank you very much for your extended response. I am not sure, will have to think about it, altough this seems right. In appreciation of your effort to help me, let me give you some color: Say I am trying to sove a 3d problem using splitting methods. Lets say that the original system f reference is xyz. One alays ends up to a system of equations in the vectorized reprezentation of the grid (very much like the array where the elements of the ma are stored). Then, when trying to solve the problem in the x direction (while in fortran storage scheme), I obtain a nice tridiagonal system of equations which I can solve very efficiently (using Thomas algorithm which is O(N) ). When I go to the second dimension, the tridiagonal system is hidden (in the original vector). However, in the rotsted yzx system it is there!! So, I am "paying the price" of writting it in the rotated view and copying it into a different vector to solve efficiently (o/wise I end up with a much more expensive sparse system to solve..) Hence the root of my questions. One of the issues is that the rows that I need to permute is rather complicated and with a little of smoke and mirrors, I hope to bypass the complications. Apologies, if all this is known to ou already or if went way too atray. Also, apologies for not catching your response earlier. Thank you agian for all your help-will look inth the references, Petros [snip] Hi Petros,
I've just uploaded: http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/li... http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/li... The *pde.cpp demonstrates 2 methods for solving pde's: 1) Alternating Directions Explicit (ADE) 2) Alternating Directions Implicit (ADI) References documenting the methods and many in source comments hopefully make it fairly easy to see how the code corresponds to the reference. If not, please let me know. The ADI method uses the solve_tridiag and rotated axes. The rotation does not rotate the axes of the array, it simply changes 1 scalar member variable in the: index_stack_length_stride_crtp_indices<> class. I hope that's what you meant by rotating axes. There are TRACE_* macros in the solve_tridiag.hpp file which illustrate what axis is actually being solve for. One thing that worries me is that the relative error for ADI is about 10x what the relative error for ADE. I've not yet figured out why :( HTH. -regards, Larry
On 05/24/11 00:58, pmamales@nyc.rr.com wrote:
Hi Larry,
Thank you very much for your extended response. I am not sure, will have to think about it, altough this seems right. In appreciation of your effort to help me, let me give you some color: Say I am trying to sove a 3d problem using splitting methods. Lets say that the original system f reference is xyz. One alays ends up to a system of equations in the vectorized reprezentation of the grid (very much like the array where the elements of the ma are stored). Then, when trying to solve the problem in the x direction (while in fortran storage scheme), I obtain a nice tridiagonal system of equations which I can solve very efficiently (using Thomas algorithm which is O(N) ). When I go to the second dimension, the tridiagonal system is hidden (in the original vector). However, in the rotsted yzx system it is there!! [snip] Hi Petros,
Based on your mention of tridiagonal system and some private emails to me, you're using the ADI method. However, Daniel Duffy, author of: http://www.amazon.com/Finite-Difference-Methods-Financial-Engineering/dp/047... expressed some doubts about ADI in this blog: http://www.datasimfinancial.com/forum/viewtopic.php?t=416 I'm a novice about PDE; so, I'd appreciate insight about why ADI seems the better solution for your problem. -regards, Larry
Hi Larry, ADI is well suited for problems with sacial and temporal vatiation of the PDE cofficients. Among experts it is indeed considered somewhat of a "hack" (not really though, only not universal solver) and other more elaborate solvers are preferred. Having said all this, it is the industry standard for its speed and ability to tackle a large variety of problems. ADE is a very new method advocated by DD and -to be honest- the fact that noone has picked it up makes me cautius (he-as of now- is the only one supporting it). Would I ever implement it? Sure, after ADI though because of the standard-ness and because once you have one you have enough infra for the other. Now on the error level you report, it does not make much sense. After all it is a 2nd order scheme, stabilized for time evolution with 2nd order RK. Are you sure that a) took the BCs into account properly? b) you propagated the BCs as well? c) have you experimented with various values of theta? d) have you experimented with different time step size? HTH P- ps: thx for the update. these days I am deep in multithreading LES (for the purpose of PDE solving) and providing with proper memory alignment for mkl to be called. -----Original Message----- From: Larry Evans Sent: Saturday, October 29, 2011 4:36 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] multi-array and pdes On 05/24/11 00:58, pmamales@nyc.rr.com wrote:
Hi Larry,
Thank you very much for your extended response. I am not sure, will have to think about it, altough this seems right. In appreciation of your effort to help me, let me give you some color: Say I am trying to sove a 3d problem using splitting methods. Lets say that the original system f reference is xyz. One alays ends up to a system of equations in the vectorized reprezentation of the grid (very much like the array where the elements of the ma are stored). Then, when trying to solve the problem in the x direction (while in fortran storage scheme), I obtain a nice tridiagonal system of equations which I can solve very efficiently (using Thomas algorithm which is O(N) ). When I go to the second dimension, the tridiagonal system is hidden (in the original vector). However, in the rotsted yzx system it is there!! [snip] Hi Petros,
Based on your mention of tridiagonal system and some private emails to me, you're using the ADI method. However, Daniel Duffy, author of: http://www.amazon.com/Finite-Difference-Methods-Financial-Engineering/dp/047... expressed some doubts about ADI in this blog: http://www.datasimfinancial.com/forum/viewtopic.php?t=416 I'm a novice about PDE; so, I'd appreciate insight about why ADI seems the better solution for your problem. -regards, Larry _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
On 10/30/11 11:28, petros wrote:
Hi Larry, ADI is well suited for problems with sacial and temporal vatiation of the PDE cofficients. Among experts it is indeed considered somewhat of a "hack" (not really though, only not universal solver) and other more elaborate solvers are preferred. Having said all this, it is the industry standard for its speed and ability to tackle a large variety of problems. ADE is a very new method advocated by DD and -to be honest- the fact that noone has picked it up makes me cautius (he-as of now- is the only one supporting it). Would I ever implement it? Sure, after ADI though because of the standard-ness and because once you have one you have enough infra for the other.
Thanks Petros for that info.
Now on the error level you report, it does not make much sense. After all it is a 2nd order scheme, stabilized for time evolution with 2nd order RK.
What is RK? The only thing I can think of is Runge-Kutta; however, I saw nothing about Runge-Kutta in Duffy's chapter 19 on which the ADI code I mentioned is based.
Are you sure that a) took the BCs into account properly? b) you propagated the BCs as well? c) have you experimented with various values of theta?
One of the theta's in Duffy's book is one of the Greeks on page 131 which is the derivative of option price w.r.t. time. The other might be the theta in equation(7.7) on p. 81; however, that's not used in chapter 19. Which theta might you be referring to?
d) have you experimented with different time step size? Yes. HTH
I've gone over the code again (and again, and again, adding comments as I go) trying to figure out what I might have done wrong. I did find that I failed to rotate the delta spaces (corrected now in the svn code. See call to std::rotate); however, I'm still getting 4-6 times as much error in innocent (the method from section 19.5 of Duffy's book) as I'm getting in Duffy's explicit. I may just try to understand the books error analysis and try to print something about that while the code is running to see why innocent is so much worse then explicit; however, that'll take some time :(
P- ps: thx for the update. these days I am deep in multithreading LES
LES = Linear Equation Solver? The innocent code uses: http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/bo... which solves a symmetric tridiagonal system. I did not use mkt or ublas or some other library because: 1) I wanted to trace the solution and that wouldn't be possible with some library. 2) After a brief look at some library(blas or something like that), it looked like it couldn't solve a set of a set of equations, it could only solve 1 set of equations (resulting in a vector result); whereas the solve_tridiag.hpp code can do that, resulting in a multi-dimensional array of results. The number of the sets of equations is the number of times the outermost loop( the one with index, i_node, in triangulate_rhs and back-substitute methods) is executed. As noted in the comments, I think this outermost loop could be made parallel; hence, multithreading could be used there.
(for the purpose of PDE solving) and providing with proper memory
As I recall, you were planning on using tridiagonal solvers also in your PDE solving. Is that still true? Are you still planning on using multi_array? If so, are you still needing to rotate the axes? If so, how do you do that with multi_array? [snip] -regards, Larry
Hi Larry, You should really try the paper, not the book (the ppaper is more recent by 4-5 years). RK is Runge-Kutta. Basically, one writes the PDE as a system of ODEs in the open set (w/out the closure). Then advancing over time is done on this system resulting in the schemes. My understanding is that the theta I am referring to is the theta of the numerical stepping, not the option sensitivity! The statement for second order time accuracy is local! According to different thetas you get the global error handling. What dimension problem are you solving? What kind of BC's. It is not clear to me whether one should solve the systems on the open or the closed set. I am of the opinion that is is the former. It should not matter, if you do Dirichlet BCs. Make sure though that you employ second order discretization at the boundary (one-sided,of course). Instead of going with the formal error analysis (always inadequate since a lot of assumptions are tacitly imposed) better do the following: Say you are solving a 2d problem. Write a function u(x,y,t) that you like. Then derive the PDE and BCs that this satisfy. Then solve with your machinery. Compare.. Yes this is nowhere near a mathematical proof but nothing in ADIs really is ;-)) (I know this is worse than others, but will halp you locate the issues: is it BCs? is it oscillations that propagate because of the one-sided derivatives? you will see). As far as milti-arrays go, I am very sceptical these days. This thing looks half=baked.. I do not trust that it can do rotated iterators. I have resorted to creating the systems to solve on the fly, by iterating on the proper plane (have created a templated functor that allows some actor class to do something in the sub-plane that I define). Typically, the LES(=linear eqtn solver) routines require copying (this is LAPAK for you) and this seems unavoidable.. The tridiagonal systems I am solving either using Thomas algorithm (wiki it), or the LAPACK tridiagonal solver, which allows for pivoting but is slower. The latter allows for many rhs vectors, which is useful only if the matrix and the BCs are independent of the dimension(s) on which you iterate. My strategy is to solve those over different threads, hence I am stuck with the plumbing. I am thinking now, if pivoting I have not looked into your code yet, partially being afraid that it would use C++ features(variadic templates) that my compiler -unfortunately- does not support. HTH, P- ps: I think this is getting overly specialized for the boost forum. Maybe we should continue privately? You have my email at home, I think. ---- Larry Evans <cppljevans@suddenlink.net> wrote:
On 10/30/11 11:28, petros wrote:
Hi Larry, ADI is well suited for problems with sacial and temporal vatiation of the PDE cofficients. Among experts it is indeed considered somewhat of a "hack" (not really though, only not universal solver) and other more elaborate solvers are preferred. Having said all this, it is the industry standard for its speed and ability to tackle a large variety of problems. ADE is a very new method advocated by DD and -to be honest- the fact that noone has picked it up makes me cautius (he-as of now- is the only one supporting it). Would I ever implement it? Sure, after ADI though because of the standard-ness and because once you have one you have enough infra for the other.
Thanks Petros for that info.
Now on the error level you report, it does not make much sense. After all it is a 2nd order scheme, stabilized for time evolution with 2nd order RK.
What is RK? The only thing I can think of is Runge-Kutta; however, I saw nothing about Runge-Kutta in Duffy's chapter 19 on which the ADI code I mentioned is based.
Are you sure that a) took the BCs into account properly? b) you propagated the BCs as well? c) have you experimented with various values of theta?
One of the theta's in Duffy's book is one of the Greeks on page 131 which is the derivative of option price w.r.t. time. The other might be the theta in equation(7.7) on p. 81; however, that's not used in chapter 19.
Which theta might you be referring to?
d) have you experimented with different time step size? Yes. HTH
I've gone over the code again (and again, and again, adding comments as I go) trying to figure out what I might have done wrong. I did find that I failed to rotate the delta spaces (corrected now in the svn code. See call to std::rotate); however, I'm still getting 4-6 times as much error in innocent (the method from section 19.5 of Duffy's book) as I'm getting in Duffy's explicit.
I may just try to understand the books error analysis and try to print something about that while the code is running to see why innocent is so much worse then explicit; however, that'll take some time :(
P- ps: thx for the update. these days I am deep in multithreading LES
LES = Linear Equation Solver? The innocent code uses:
http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/bo...
which solves a symmetric tridiagonal system. I did not use mkt or ublas or some other library because:
1) I wanted to trace the solution and that wouldn't be possible with some library.
2) After a brief look at some library(blas or something like that), it looked like it couldn't solve a set of a set of equations, it could only solve 1 set of equations (resulting in a vector result); whereas the solve_tridiag.hpp code can do that, resulting in a multi-dimensional array of results. The number of the sets of equations is the number of times the outermost loop( the one with index, i_node, in triangulate_rhs and back-substitute methods) is executed. As noted in the comments, I think this outermost loop could be made parallel; hence, multithreading could be used there.
(for the purpose of PDE solving) and providing with proper memory
As I recall, you were planning on using tridiagonal solvers also in your PDE solving. Is that still true? Are you still planning on using multi_array? If so, are you still needing to rotate the axes? If so, how do you do that with multi_array? [snip]
-regards, Larry
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
participants (4)
-
Larry Evans -
petros -
Pierre-Andre Noel -
pmamales＠nyc.rr.com