Contents

Step 1: Precondition and Postcondition
Step 2: Determine Loop-Invariants
Step 3: Determine the Loop-Guard
Step 4: Determine the Initialization
Step 5: Progressing through the Matrices
Step 6: State after Repartitioning
Step 7: State after Moving the Lines
Step 8: Determining the Update
Final Algorithm

Derivation of LU Factorization Unblocked Variant 1

Operation to be computed:

$\alpha \becomes APDOT(x, y, \alpha)$

where $x$ and $y$ are vectors of the same length, and $\alpha$ is their dot product (a scalar). We start with the blank worksheet

- $\resetsteps \footnotesize \worksheet$

Step 1: Precondition and Postcondition

The precondition is given by the predicate

$\alpha = \hat{\alpha},$

where $\hat{\alpha}$ denotes the original contents of $ A $ and is a dummy variable. The postcondition is given by the predicate

$\alpha = x^T y + \hat{\alpha}.$

These conditions are entered in the worksheet as Step 1a and 1b, as in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \color{red} \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \color{red} \alpha = \tr{x} y + \hat{\alpha} } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 2: Determine Loop-Invariants

As the computation progresses, the algorithms will sweep through the vectors exposing partitions:

$\begin{eqnarray*} x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} &,& y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} , \\ \end{eqnarray*}$

where $x_{L}$ and $y_{L}$ are both vectors and of equal size.

To determine loop-invariants, we take these partitioned vectors and substitute them into the postcondition:

$\begin{eqnarray*} \alpha &=& \FlaOneByTwo{x_{T}}{x_{B}} \FlaTwoByOne{y_{T}} {y_{B}} ~ + ~ \hat{\alpha} \end{eqnarray*}$

Manipulating this term yeilds

$\alpha = x_{T}^T y_{T} + y_{B}^T y_{B} + \hat{\alpha}.$

Depending on whether we choose to start from the top or the bottom of the vectors $x$ and $y$ , we find that $\alpha$ can, at an intermediate stage, be in one of the following two states:

Invariant 1	$\alpha = \tr{x}_T y_T + \hat{\alpha}$
Invariant 2	$\alpha = \tr{x}_B y_B + \hat{\alpha}$

We will pick Invariant 1 for our subsequent discussion, which is entered in the worksheet as in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ \color{red} ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 3: Determine the Loop-Guard

Next we determine the loop-guard $ G $ . Note that

$\left( \alpha = \tr{x}_T y_T + \hat{\alpha} ~ \wedge ~ 0 \leq m( x_T ) = m( y_T ) \leq m( x ) \right) \wedge \neg G$

must imply the postcondition. Taking $G = ( m( x_{T} ) < m( x ) )$ has this property.

The loop-guard is entered as Step 3, in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } % Step 3: Loop-guard \renewcommand{\guard}{ \color{red} m( x_T ) < m( x ) } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 4: Determine the Initialization

Now we determine an initialization that puts $\alpha$ in a state where the loop-invariant is true. The partitioning

$\begin{eqnarray*} x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} &,& y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} , \\ \end{eqnarray*}$

with $x_{T}$ and $\hat{y}_{T}$ both empty ( $0 \times 0$ ) has the desired property. This is entered as Step 4, in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } % Step 3: Loop-guard \renewcommand{\guard}{ m( x_T ) < m( x ) } % Step 4: Initialize \renewcommand{\partitionings}{ \color{red} $ x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} $ , $ y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} $ } \renewcommand{\partitionsizes}{ \color{red} $ x_{T} $ and $ y_{T} $ have $ 0 $ elements } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 5: Progressing through the Matrices

The initialization sets the top partitions to $0 \times 0$ while the loop-guard indicates that the loop completes when these same paritions envelop the entire vectors. Thus, the algorithm must march through the vectors in a way that grows the top partition. Since we are trying to derive an unblocked algorithm, we choose to expose columns at the top of the loop that are then moved to the other side of the thick lines at the bottom of the loop, as in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } % Step 3: Loop-guard \renewcommand{\guard}{ m( x_T ) < m( x ) } % Step 4: Initialize \renewcommand{\partitionings}{ $ x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} $ , $ y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} $ } \renewcommand{\partitionsizes}{ $ x_{T} $ and $ y_{T} $ have $ 0 $ elements } % Step 5a: Repartition the operands \renewcommand{\repartitionings}{ \color{red} $ \FlaTwoByOne{ x_T } { x_B } \rightarrow \FlaThreeByOneB{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \rightarrow \FlaThreeByOneB{y_0} {\psi_1} {y_2} $ } \renewcommand{\repartitionsizes}{ \color{red} $ \chi_1 $ and $ \psi_1 $ are scalars} % Step 5b: Move the double lines \renewcommand{\moveboundaries}{ \color{red} $ \FlaTwoByOne{ x_T } { x_B } \leftarrow \FlaThreeByOneT{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \leftarrow \FlaThreeByOneT{y_0} {\psi_1} {y_2} $ } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 6: State after Repartitioning

The statement Repartition ... with ... exposes subvectors. We determine what the state of $\alpha$ is in terms of these subvectors via textual substitution into the loop-invariant.

The repartitioning step

$\begin{eqnarray*} \FlaTwoByOne{ x_T } { x_B } \rightarrow \FlaThreeByOneB{x_0} {\chi_1} {x_2} &,& \FlaTwoByOne{ y_T } { y_B } \rightarrow \FlaThreeByOneB{y_0} {\psi_1} {y_2} , \\ \end{eqnarray*}$

can be substitued into the loop-invariant, yielding

$\begin{eqnarray*} \alpha &=& x_0^T y_0 + \hat{\alpha} ~ \wedge ~ 0 \leq m( x_0 ) = m( y_0 ) \leq m( x ) \\ \end{eqnarray*}$

This is then entered in Step 6, in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } % Step 3: Loop-guard \renewcommand{\guard}{ m( x_T ) < m( x ) } % Step 4: Initialize \renewcommand{\partitionings}{ $ x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} $ , $ y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} $ } \renewcommand{\partitionsizes}{ $ x_{T} $ and $ y_{T} $ have $ 0 $ elements } % Step 5a: Repartition the operands \renewcommand{\repartitionings}{ $ \FlaTwoByOne{ x_T } { x_B } \rightarrow \FlaThreeByOneB{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \rightarrow \FlaThreeByOneB{y_0} {\psi_1} {y_2} $ } \renewcommand{\repartitionsizes}{ $ \chi_1 $ and $ \psi_1 $ are scalars} % Step 5b: Move the double lines \renewcommand{\moveboundaries}{ $ \FlaTwoByOne{ x_T } { x_B } \leftarrow \FlaThreeByOneT{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \leftarrow \FlaThreeByOneT{y_0} {\psi_1} {y_2} $ } % Step 6: State after repartitioning % Note: The below needs editing!!! \renewcommand{\beforeupdate}{ {\color{red} ( \alpha = \tr{x}_0 y_0 + \hat{\alpha} ) \wedge ( 0 \leq m( x_0 ) = m( y_0 ) < m( x ) ) } } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 7: State after Moving the Lines

The statement Continue with ... redefines the quadrants. We determine what the state of matrix must be if the loop-invariant is to again hold at the bottom of the loop, again via textual substitution into the loop-invariant.

The Continue with ... step The repartitioning step

$\begin{eqnarray*} \FlaTwoByOne{ x_T } { x_B } \leftarrow \FlaThreeByOneT{x_0} {\chi_1} {x_2} &,& \FlaTwoByOne{ y_T } { y_B } \leftarrow \FlaThreeByOneT{y_0} {\psi_1} {y_2} , \\ \end{eqnarray*}$

Substituting these into the loop-invariant yields

$\begin{eqnarray*} \alpha &=& x_0^T y_0 + \chi_1 \psi_1 + \hat{\alpha}. \\ \end{eqnarray*}$

This is then entered in Step 7, in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } % Step 3: Loop-guard \renewcommand{\guard}{ m( x_T ) < m( x ) } % Step 4: Initialize \renewcommand{\partitionings}{ $ x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} $ , $ y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} $ } \renewcommand{\partitionsizes}{ $ x_{T} $ and $ y_{T} $ have $ 0 $ elements } % Step 5a: Repartition the operands \renewcommand{\repartitionings}{ $ \FlaTwoByOne{ x_T } { x_B } \rightarrow \FlaThreeByOneB{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \rightarrow \FlaThreeByOneB{y_0} {\psi_1} {y_2} $ } \renewcommand{\repartitionsizes}{ $ \chi_1 $ and $ \psi_1 $ are scalars} % Step 5b: Move the double lines \renewcommand{\moveboundaries}{ $ \FlaTwoByOne{ x_T } { x_B } \leftarrow \FlaThreeByOneT{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \leftarrow \FlaThreeByOneT{y_0} {\psi_1} {y_2} $ } % Step 6: State after repartitioning % Note: The below needs editing!!! \renewcommand{\beforeupdate}{ ( \alpha = \tr{x}_0 y_0 + \hat{\alpha} ) \wedge ( 0 \leq m( x_0 ) = m( y_0 ) < m( x ) ) } % Step 7: State after moving of double lines % Note: The below needs editing!!! \renewcommand{\afterupdate}{ { \color{red} ( \alpha = \tr{x}_0 y_0 + \chi_1 \psi_1 + \hat{\alpha} ) \wedge ( 0 < m( x_0 ) + 1 = m( y_0 ) + 1 \leq m( x ) ) } } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Step 8: Determining the Update

The question now becomes how to update the contents of matrix $\alpha$ from

$x_0^T y_0 + \hat{\alpha}$

$x_0^T y_0 + \chi_1 \psi_1 + \hat{\alpha}$

This justifies the update

$\begin{array}{l} \alpha \becomes \chi_1 \psi_1 + \alpha \end{array}$

which is entered as Step 8, in

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\operation}{ \alpha \becomes \tr{x} y + \alpha } \renewcommand{\routinename}{ \alpha \becomes \mbox{\sc sapdot}(x, y, \alpha) } \renewcommand{\routinecost}{ 2m } % Step 1a: Precondition \renewcommand{\precondition}{ \alpha = \hat{\alpha} \wedge 0 \leq m( x ) = m( y ) } % Step 1b: Postcondition \renewcommand{\postcondition}{ \alpha = \tr{x} y + \hat{\alpha} } % Step 2: Invariant % Note: Right-hand side of equalities must be updated appropriately \renewcommand{\invariant}{ ( \alpha = \tr{x}_T y_T + \hat{\alpha} ) \wedge ( 0 \leq m( x_T ) = m( y_T ) \leq m( x ) ) } % Step 3: Loop-guard \renewcommand{\guard}{ m( x_T ) < m( x ) } % Step 4: Initialize \renewcommand{\partitionings}{ $ x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} $ , $ y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} $ } \renewcommand{\partitionsizes}{ $ x_{T} $ and $ y_{T} $ have $ 0 $ elements } % Step 5a: Repartition the operands \renewcommand{\repartitionings}{ $ \FlaTwoByOne{ x_T } { x_B } \rightarrow \FlaThreeByOneB{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \rightarrow \FlaThreeByOneB{y_0} {\psi_1} {y_2} $ } \renewcommand{\repartitionsizes}{ $ \chi_1 $ and $ \psi_1 $ are scalars} % Step 5b: Move the double lines \renewcommand{\moveboundaries}{ $ \FlaTwoByOne{ x_T } { x_B } \leftarrow \FlaThreeByOneT{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \leftarrow \FlaThreeByOneT{y_0} {\psi_1} {y_2} $ } % Step 6: State after repartitioning % Note: The below needs editing!!! \renewcommand{\beforeupdate}{ ( \alpha = \tr{x}_0 y_0 + \hat{\alpha} ) \wedge ( 0 \leq m( x_0 ) = m( y_0 ) < m( x ) ) } % Step 7: State after moving of double lines % Note: The below needs editing!!! \renewcommand{\afterupdate}{ ( \alpha = \tr{x}_0 y_0 + \chi_1 \psi_1 + \hat{\alpha} ) \wedge ( 0 < m( x_0 ) + 1 = m( y_0 ) + 1 \leq m( x ) ) } % Step 8: Insert the updates required to change the % state from that given in Step 6 to that given in Step 7 % Note: The below needs editing!!! \renewcommand{\update}{ \color{red} $ % \begin{tabular}{l @{\hspace{10ex}} | @{\hspace{10ex}} l} % {\bf Update:} & {\bf Cost:} \\ \alpha \becomes \chi_1 \psi_1 + \alpha % & ~~~~~~$2$ flops % \end{tabular} $ } \worksheet % this command generates the worksheet using the % commands that were defined between the \resetsteps % and the worksheet command$

Final Algorithm

By eliminating all references to the dummy variables $\hat{\alpha}$ , the algorithm for computing Apdot is left:

$\resetsteps % Reset all the commands to create a blank worksheet % Define the operation to be computed \renewcommand{\routinename}{ \alpha \becomes \hat{\alpha} + x^T y } % Step 3: Loop-guard \renewcommand{\guard}{ m( x_T ) < m( x ) } % Step 4: Initialize \renewcommand{\partitionings}{ $ x \rightarrow \FlaTwoByOne{x_{T}} {x_{B}} $ , $ y \rightarrow \FlaTwoByOne{y_{T}} {y_{B}} $ } \renewcommand{\partitionsizes}{ $ x_{T} $ and $ y_{T} $ have $ 0 $ elements } % Step 5a: Repartition the operands \renewcommand{\repartitionings}{ $ \FlaTwoByOne{ x_T } { x_B } \rightarrow \FlaThreeByOneB{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \rightarrow \FlaThreeByOneB{y_0} {\psi_1} {y_2} $ } \renewcommand{\repartitionsizes}{ $ \chi_1 $ and $ \psi_1 $ are scalars} % Step 5b: Move the double lines \renewcommand{\moveboundaries}{ $ \FlaTwoByOne{ x_T } { x_B } \leftarrow \FlaThreeByOneT{x_0} {\chi_1} {x_2} $ , $ \FlaTwoByOne{ y_T } { y_B } \leftarrow \FlaThreeByOneT{y_0} {\psi_1} {y_2} $ } % Step 8: Insert the updates required to change the % state from that given in Step 6 to that given in Step 7 % Note: The below needs editing!!! \renewcommand{\update}{ $ \begin{array}{l} \alpha \becomes \chi_1 \psi_1 + \alpha \end{array} $ } \footnotesize \FlaAlgorithmNarrow % this command generates the algorithm$

LinearAlgebraWiki: Apdot/Derivation/UnbVar1 (last edited 2007-05-24 22:13:51 by MarthaGanser)