blitz_3.html

c++经典教材 Blitz++ v0.8

<TABLE><tr><td>&nbsp;</td><td class=smallexample><FONT SIZE=-1><pre>#include &#60;blitz/array.h&#62;

using namespace blitz;

int main()
{
    Array&#60;float,1&#62; x(4), y(4);
    Array&#60;float,2&#62; A(4,4);

    x = 1, 2, 3, 4;
    y = 1, 0, 0, 1;

    firstIndex i;
    secondIndex j;

    A = x(i) * y(j);

    cout &#60;&#60; A &#60;&#60; endl;

    return 0;
}

</FONT></pre></td></tr></table></P><P>

And the output:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=smallexample><FONT SIZE=-1><pre>4 x 4
[         1         0         0         1 
          2         0         0         2 
          3         0         0         3 
          4         0         0         4 ]

</FONT></pre></td></tr></table></P><P>

Index placeholders can <EM>not</EM> be used on the left-hand side of an
expression.  If you need to reorder the indices, you must do this on the
right-hand side.
</P><P>

In real-world tensor notation, repeated indices imply a contraction (or
summation).  For example, this tensor expression computes a matrix-matrix
product:
<pre>
 ij    ik  kj
C   = A   B
</pre>
</P><P>

The repeated k index is interpreted as meaning
<pre>
c    = sum of {a   * b  } over k
 ij             ik    kj
</pre>
</P><P>

<A NAME="IDX285"></A>
<A NAME="IDX286"></A>
</P><P>

In Blitz++, repeated indices do <EM>not</EM> imply contraction.  If you want
to contract (sum along) an index, you must use the <CODE>sum()</CODE> function:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,2&#62; A, B, C;   // ...
firstIndex i;
secondIndex j;
thirdIndex k;

C = sum(A(i,k) * B(k,j), k);
</pre></td></tr></table></P><P>

The <CODE>sum()</CODE> function is an example of an <EM>array reduction</EM>,
described in the next section.
</P><P>

Index placeholders can be used in any order in an expression.  This example
computes a kronecker product of a pair of two-dimensional arrays, and
permutes the indices along the way:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,2&#62; A, B;   // ...
Array&#60;float,4&#62; C;      // ...
fourthIndex l;

C = A(l,j) * B(k,i);
</pre></td></tr></table></P><P>

This is equivalent to the tensor notation
<pre>
 ijkl    lj ki
C     = A  B
 </pre>
</P><P>

Tensor-like notation can be mixed with other array notations:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,2&#62; A, B;  // ...
Array&#60;double,4&#62; C;    // ...

C = cos(A(l,j)) * sin(B(k,i)) + 1./(i+j+k+l);
</pre></td></tr></table></P><P>

<A NAME="IDX287"></A>
An important efficiency note about tensor-like notation: the right-hand side
of an expression is <EM>completely evaluated</EM> for <EM>every</EM> element in
the destination array.  For example, in this code:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,1&#62; x(4), y(4);
Array&#60;float,2&#62; A(4,4):

A = cos(x(i)) * sin(y(j));
</pre></td></tr></table></P><P>

The resulting implementation will look something like this:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>for (int n=0; n &#60; 4; ++n)
  for (int m=0; m &#60; 4; ++m)
    A(n,m) = cos(x(n)) * sin(y(m));
</pre></td></tr></table></P><P>

The functions <CODE>cos</CODE> and <CODE>sin</CODE> will be invoked sixteen times each.
It's possible that a good optimizing compiler could hoist the <CODE>cos</CODE>
evaluation out of the inner loop, but don't hold your breath -- there's a
lot of complicated machinery behind the scenes to handle tensor notation,
and most optimizing compilers are easily confused.  In a situation like the
above, you are probably best off manually creating temporaries for
<CODE>cos(x)</CODE> and <CODE>sin(y)</CODE> first.
</P><P>

<HR SIZE="6">
<A NAME="SEC94"></A>
<TABLE CELLPADDING=1 CELLSPACING=1 BORDER=0>
<TR><TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC93"> &lt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC95"> &gt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC76"> &lt;&lt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC74"> Up </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_4.html#SEC98"> &gt;&gt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz.html#SEC_Top">Top</A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_toc.html#SEC_Contents">Contents</A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[Index]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_abt.html#SEC_About"> ? </A>]</TD>
</TR></TABLE>
<H2> 3.12 Array reductions </H2>
<!--docid::SEC94::-->
<P>

Currently, Blitz++ arrays support two forms of reduction:
</P><P>

<UL>

<LI>Reductions which transform an array into a scalar (for example,
summing the elements).  These are referred to as <STRONG>complete
reductions</STRONG>.
<P>

<LI>Reducing an N dimensional array (or array expression) to an N-1
dimensional array expression.  These are called <STRONG>partial reductions</STRONG>.
<P>

</UL>
<P>

<A NAME="IDX288"></A>
<A NAME="IDX289"></A>
<A NAME="IDX290"></A>
</P><P>

<HR SIZE="6">
<A NAME="SEC95"></A>
<TABLE CELLPADDING=1 CELLSPACING=1 BORDER=0>
<TR><TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC94"> &lt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC96"> &gt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC76"> &lt;&lt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC74"> Up </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_4.html#SEC98"> &gt;&gt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz.html#SEC_Top">Top</A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_toc.html#SEC_Contents">Contents</A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[Index]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_abt.html#SEC_About"> ? </A>]</TD>
</TR></TABLE>
<H2> 3.13 Complete reductions </H2>
<!--docid::SEC95::-->
<P>

Complete reductions transform an array (or array expression) into 
a scalar.  Here are some examples:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,2&#62; A(3,3);
A = 0, 1, 2,
    3, 4, 5,
    6, 7, 8;
cout &#60;&#60; sum(A) &#60;&#60; endl          // 36
     &#60;&#60; min(A) &#60;&#60; endl          // 0
     &#60;&#60; count(A &#62;= 4) &#60;&#60; endl;  // 5
</pre></td></tr></table></P><P>

Here are the available complete reductions:
</P><P>

<DL COMPACT>
<DT><CODE>sum()</CODE>
<DD><A NAME="IDX291"></A>
Summation (may be promoted to a higher-precision type)
<P>

<DT><CODE>product()</CODE>
<DD><A NAME="IDX292"></A>
Product 
<P>

<DT><CODE>mean()</CODE>
<DD><A NAME="IDX293"></A>
Arithmetic mean (promoted to floating-point type if necessary) 
<P>

<DT><CODE>min()</CODE>
<DD><A NAME="IDX294"></A>
Minimum value 
<P>

<DT><CODE>max()</CODE>
<DD><A NAME="IDX295"></A>
Maximum value 
<P>

<DT><CODE>minIndex()</CODE>
<DD><A NAME="IDX296"></A>
Index of the minimum value (<CODE>TinyVector&#60;int,N_rank&#62;</CODE>)
<P>

<DT><CODE>maxIndex()</CODE>
<DD><A NAME="IDX297"></A>
Index of the maximum value (<CODE>TinyVector&#60;int,N_rank&#62;</CODE>)
<P>

<DT><CODE>count()</CODE>
<DD><A NAME="IDX298"></A>
Counts the number of times the expression is logical true (<CODE>int</CODE>)
<P>

<DT><CODE>any()</CODE>
<DD><A NAME="IDX299"></A>
True if the expression is true anywhere (<CODE>bool</CODE>)
<P>

<DT><CODE>all()</CODE>
<DD><A NAME="IDX300"></A>
True if the expression is true everywhere (<CODE>bool</CODE>)
</DL>
<P>

<STRONG>Note:</STRONG> <CODE>minIndex()</CODE> and <CODE>maxIndex()</CODE> return TinyVectors, 
even when the rank of the array (or array expression) is 1.
</P><P>

Reductions can be combined with <CODE>where</CODE> expressions (<A HREF="blitz_3.html#SEC97">3.15 where statements</A>)
to reduce over some part of an array.  For example, <CODE>sum(where(A &#62; 0,
A, 0))</CODE> sums only the positive elements in an array.
</P><P>

<HR SIZE="6">
<A NAME="SEC96"></A>
<TABLE CELLPADDING=1 CELLSPACING=1 BORDER=0>
<TR><TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC95"> &lt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC97"> &gt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC76"> &lt;&lt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_3.html#SEC74"> Up </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_4.html#SEC98"> &gt;&gt; </A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT"> &nbsp; <TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz.html#SEC_Top">Top</A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_toc.html#SEC_Contents">Contents</A>]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[Index]</TD>
<TD VALIGN="MIDDLE" ALIGN="LEFT">[<A HREF="blitz_abt.html#SEC_About"> ? </A>]</TD>
</TR></TABLE>
<H2> 3.14 Partial Reductions </H2>
<!--docid::SEC96::-->
<P>

<A NAME="IDX301"></A>
<A NAME="IDX302"></A>
<A NAME="IDX303"></A>
</P><P>

Here's an example which computes the sum of each row of a two-dimensional
array:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,2&#62; A;    // ...
Array&#60;float,1&#62; rs;   // ...
firstIndex i;
secondIndex j;

rs = sum(A, j);
</pre></td></tr></table></P><P>

The reduction <CODE>sum()</CODE> takes two arguments:
</P><P>

<UL>

<LI>The first argument is an array or array expression.
<P>

<LI>The second argument is an index placeholder indicating the
dimension over which the reduction is to occur.  
<P>

</UL>
<P>

Reductions have an <STRONG>important restriction</STRONG>: It is currently only
possible to reduce over the <EM>last</EM> dimension of an array or array
expression.  Reducing a dimension other than the last would require Blitz++
to reorder the dimensions to fill the hole left behind.  For example, in
order for this reduction to work:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>Array&#60;float,3&#62; A;   // ...
Array&#60;float,2&#62; B;   // ...
secondIndex j;

// Reduce over dimension 2 of a 3-D array?
B = sum(A, j);
</pre></td></tr></table></P><P>

Blitz++ would have to remap the dimensions so that the third dimension
became the second.  It's not currently smart enough to do this.
</P><P>

However, there is a simple workaround which solves some of the problems
created by this limitation: you can do the reordering manually, prior to the
reduction:
</P><P>

<TABLE><tr><td>&nbsp;</td><td class=example><pre>B = sum(A(i,k,j), k);
</pre></td></tr></table></P><P>

Writing <CODE>A(i,k,j)</CODE> interchanges the second and third dimensions,
permitting you to reduce over the second dimension.  Here's a list of the
reduction operations currently supported:
</P><P>

<DL COMPACT>
<DT><CODE>sum()</CODE>
<DD>Summation
<P>

<DT><CODE>product()</CODE>
<DD>Product 
<P>

<DT><CODE>mean()</CODE>
<DD>Arithmetic mean (promoted to floating-point type if necessary)
<P>

<DT><CODE>min()</CODE>
<DD>Minimum value
<P>

<DT><CODE>max()</CODE>
<DD>Maximum value
<P>

<DT><CO