mkl_?bsrmm documentation question

[Nicholas_K_2]

mkl_?bsrmm is:C := alpha op(A) * B + beta * C, in bsr format)

For the dense matrix B (often tall and skinny, eg block Krylov methods), the documentation specifies:

b: [...types...] Array,DIMENSION(ldb,n)for one-based indexing,DIMENSION(m,ldb)for zero-based indexing.

On entry withtransa='N'or'n', the leadingn-by-kblock part of the arraybmust contain the matrixB, otherwise the leadingm-by-nblock part of the arraybmust contain the matrixB.

ldb:INTEGER. Specifies the first dimension (in blocks) ofbas declared in the calling (sub)program.

A few questions.My use case: transa='n', beta=0, zero-based indexing.

1) Since A is m-by-k and C is m-by-n then B must bek-by-n... right?

2) Assuming this, what does it mean for b to have DIMENSION(m, ldb)? Often in my program m < k, and I cannot satisfy this constraint. Yet, the documentation does not explicitly exclude "short fat" matrices (m < k).

3) Is ldb the first or second dimension of b?

4)In my use case, B and Care submatrices of larger, row-major matrices. Do I need to repack in order to use this routine?

Thanks for your help!

--Nick

[Chao_Y_Intel · Accepted solution]

1>Since A is m-by-k and C is m-by-n then B must be k-by-n ... right?

Yes, It looks the document error in the manual

4>B and C are submatrices of larger, row-major matrices. Do I need to repack in order to use this routine?

The functions support row-major function. Check the matdescra parameter setting, and some description on the parameter:
http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-34C8DB79-0139-46E0-8B53-99F3BEE7B2D4.htm#TBL2-6

You needs set it as zero-based indexing

2) Assuming this, what does it mean for b to have DIMENSION(m, ldb)? Often in my program m < k, and I cannot satisfy this constraint. Yet, the documentation does not explicitly exclude "short fat" matrices (m < k).
For zero-based, it is (k, ldb) actually.

3) Is ldb the first or second dimension of b?
For zero-based, it is the first dimension.

Thanks,
Chao

[Nicholas_K_2]

Chao - thanks for your reply! I'm excited to hear that this routine will work for me.

Two clarifications and two more questions:

1)In the matdescra documentation you linked, the only line I can find that supports your statement is

In Fortran the column-major order is used, and in C - row-major order.

Does this meanI must use the C Sparse BLAS interface in order for MKL to interpret my dense arrays as row-major, or canIalso use the Fortran interface with 0-based indexing?

2) You said:

3) Is ldb the first or second dimension of b?

For zero-based, it is the first dimension.

This means that ldb is the leading dimension of b, and since b is stored row-major, ldb gives the number of columns in b, correct?

3) Are all the dimension parameters (m,n,k,ldb,ldc) in units of the block size lb? The documentation says

n:INTEGER. Number of columns of the matrixC.

but the other parameters (m,k,ldb,ldc) all specify block rows/columns.

4) In the BSR format with 0-based indexing, the indx and val arrays are stored (sparse) block-row-major with (dense) row-major blocks. Does this mean that b and cmust be stored (dense) block-row-major with (dense) row-major blocks? Or can they be stored in (unblocked, dense) row-major layout, but zero-padded to a multiple of lb?

Thanks again for your help.

--Nick

[Sergey_K_Intel1]

Dear Nick,

1. Does this meanI must use the C Sparse BLAS interface in order for MKL to interpret my dense arrays as row-major, or canIalso use the Fortran interface with 0-based indexing?

Ifb and c in your Fortran program are stored in row-major fashion,you can also use the 0-based switch from Fortran.

2. This means that ldb is the leading dimension of b, and since b is stored row-major, ldb gives the number of columns in b, correct?

More precisely ldb specifies the leading dimension of b for one-based indexing, and the second dimension of b for zero-based indexing, as declared in the calling (sub)program.ldb must be >= m*lb for the non-transposed case and >=kb*lb otherwise in the case of 1-based indexing. If we choose 0-based indexing ldb >= k*lb for non-transposed A, and ldb >= m*lb for transposed.

3.Are all the dimension parameters (m,n,k,ldb,ldc) in units of the block size lb?

m and k are onlyblock dimensions. The rest of parameters are just unblocked dimensions.


4) In the BSR format with 0-based indexing, the indx and val arrays are stored (sparse) block-row-major with (dense) row-major blocks. Does this mean that b and cmust be stored (dense) block-row-major with (dense) row-major blocks? Or can they be stored in (unblocked, dense) row-major layout, but zero-padded to a multiple of lb?

b and c must be stored in regular unblocked denseway as for the rest MKL Sparse ordenseBLAS . If youchoose 1-based indexing, B and Care strored in column-major fashion and row-major otherwise.

As concerns matrix A, blocks are always listed in row-major fashion but elements in blocks are stored in row-major layout for 0-based indexing and column-major for 1-based indexing. Please see Appendix Ain MKL Reference Manual for more detailsabout the BSR storage.

All the best
Sergey

[Nicholas_K_2]

Sergey says:

"More precisely ldb specifies the leading dimension of b for one-based indexing, and the second dimension of b for zero-based indexing, as declared in the calling (sub)program.ldb must be >= m*lb for the non-transposed case and >=kb*lb otherwise in the case of 1-based indexing. If we choose 0-based indexing ldb >= k*lb for non-transposed A, and ldb >= m*lb for transposed."

I understand what you are saying for the case of 0-based indexing, but not for 1-based indexing (bold):

• When performing A*B, when A is m-by-k, why would either dimension of B be constrained bym?
• Likewise, when performing A^T*B, why would either parameter of B depend on k?

It seems to me that it would be

• transa='n' ==> ldb >= k*lb
• transa='y' ==> ldb >= m*lb

for both 0-based and 1-based indexing... Sorry for my confusion; thanks for your help! Best, --Nick

[mecej4]

When a transpose is to be used, that is taken into account after entering the MKL routine(s).

Since the declared sizes of the array arguments passed may be larger than the actual sizes, the information concerning the declared size needs to available to MKL regardless of whether the transpose is specified for one or more arguments.

In order to compute the effective 1-D array index for a 2-D array, in Fortran one needs the leading dimension, i.e., the first dimension in the declaration of the array. In C, one needs the trailing dimension.

Whether the arrays are zero-based or one-based, or even start with a negative or positive base, is merely a distraction. After all, matrix operations are defined in mathematics with the convention of rows and columns being numbered starting from 1. Think of the additional argument as indicating the number of rows or columns, rather than the sequence number of the last row or column. The former is invariant with respect to changing from 1-base to 0-base numbering; the latter is not.