I have some Fortran code which contains calls to Lapack and Blas functions that I am trying to compile using f2py with:
- Windows 10
- Numpy version 1.19
- Anaconda Python 3.8.5
The clean version of the code I am trying to compile can be found in the link: https://software.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook -samples-120115.zip under /BlockTDS_SPD/source/.
An explanation of this function can be found in this link: https://software.intel.com/content/www/us/en/develop/documentation/mkl-cookbook/top/factoring-block-tridiagonal-symmetric-positive-definite-matrices.html .
Here is the code that I am trying to compile with Cf2py intent statements added:
!***********************************************************************
! Copyright(C) 2014-2015 Intel Corporation. All Rights Reserved.
!
! The source code, information and material ("Material") contained herein is
! owned by Intel Corporation or its suppliers or licensors, and title to such
! Material remains with Intel Corporation or its suppliers or licensors. The
! Material contains proprietary information of Intel or its suppliers and
! licensors. The Material is protected by worldwide copyright laws and treaty
! provisions. No part of the Material may be used, copied, reproduced,
! modified, published, uploaded, posted, transmitted, distributed or disclosed
! in any way without Intel's prior express written permission. No license
! under any patent, copyright or other intellectual property rights in the
! Material is granted to or conferred upon you, either expressly, by
! implication, inducement, estoppel or otherwise. Any license under such
! intellectual property rights must be express and approved by Intel in
! writing.
!
! *Third Party trademarks are the property of their respective owners.
!
! Unless otherwise agreed by Intel in writing, you may not remove or alter
! this notice or any other notice embedded in Materials by Intel or Intel's
! suppliers or licensors in any way.
!
!***********************************************************************
! Content:
! Subroutine DPBLTRF for Cholesky factorization of symmetric
! positive definite block tridiagonal matrix.
!***********************************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Definition:
! ===========
! SUBROUTINE DPBLTRF(N, NB, D, LDD, B, LDB, INFO)
!
! ..Scalar arguments..
! INTEGER N, NB, LDD, LDB, INFO
! ..
! ..Array arguments..
! REAL*8 D(LDD,*), B(LDB,*)
! ..
! Purpose:
! ========
! DPBLTRF computes Cholesky L*L^t-factorization of symmetric positive
! definite block tridiagonal matrix A
! D_1 B_1^t
! B_1 D_2 B_2^t
! B_2 D_3 B_3^t
! . . .
! . . .
! B_N-2 D_N-1 B_N-1^t
! B_N-1 D_N
! The factorization has the form A = L*L**t, where L is a lower
! bidiagonal block matrix
! L_1
! C_1 L_2
! C_2 L_3
! . . .
! . . .
! C_N-2 L_N-1
! C_N-1 L_N
! This is a block version of LAPACK DPTTRF subroutine.
!
! Arguments:
! ==========
! N (input) INTEGER
! The number of block rows of the matrix A. N >= 0.
!
! NB (input) INTEGER
! The size of blocks. NB >= 0.
!
! D (input/output) REAL*8 array, dimension (LDD,N*NB)
! On entry, the array stores N diagonal blocks (each of size NB by
! NB) of the matrix to be factored. The blocks are stored
! sequentially: first NB columns of D store block D_1, second NB
! columns store block D_2,...,last NB columns store block D_N.
! Note: As the diagonal blocks are symmetric only lower or upper
! ====
! triangle is needed to store blocks' elements. In this code
! lower storage is used!!!
! On exit, the array stores diagonal blocks of triangular factor L.
! Diagonal blocks of lower triangular factor L replace
! respective lower triangles of blocks D_j (1 <= j <= N).
! Caution: upper triangles of diagonal blocks are not zeroed on
! =======
! exit!!!
!
! LDD (input) INTEGER
! The leading dimension of array D. LDD >= NB.
!
! B (input/output) REAL*8 array, dimension (LDB,(N-1)*NB)
! On entry, the array stores sub-diagonal blocks (each of size NB
! by NB) of the matrix to be factored. The blocks are stored
! sequentially: first NB columns of B store block B_1, second
! NB columns store block B_2,...,last NB columns store block
! B_N-1.
! On exit, the array stores sub-diagonal blocks of triangular factor
! L.
!
! LDB (input) INTEGER
! The leading dimension of array B. LDB >= NB.
!
! INFO (output) INTEGER
! = 0: successful exit
! < 0: if INFO = -i, the i-th argument had an illegal value
! > 0: if INFO = i, the leading minor of order i (and
! therefore the matrix A itself) is not
! positive-definite, and the factorization could not be
! completed. This may indicate an error in forming the
! matrix A.
! =====================================================================
SUBROUTINE DPBLTRF(N, NB, D, LDD, B, LDB, INFO)
IMPLICIT NONE
! ..Scalar arguments..
INTEGER N, NB, LDD, LDB, INFO
! ..Array arguments..
REAL*8 D(LDD,*), B(LDB,*)
! =====================================================================
! .. Local Scalars ..
INTEGER K
Cf2py integer, intent(in) N, NB, LDD, LDB
Cf2py intent(in,out,copy) :: D, B
Cf2py integer, intent(out) :: INFO
! ..
! .. Executable Statements ..
! ..
! Test the input arguments.
INFO = 0
IF(N .LT. 0) THEN
INFO = -1
ELSE IF(NB .LT. 0) THEN
INFO = -2
ELSE IF(LDD .LT. NB) THEN
INFO = -4
ELSE IF(LDB .LT. NB) THEN
INFO = -6
END IF
IF(INFO .NE. 0) THEN
RETURN
END IF
! ..
! Compute Cholesky factorization of the first diagonal block
CALL DPOTRF('L', NB, D, LDD, INFO)
IF(INFO .NE. 0) THEN
RETURN
END IF
!
! Main loop
DO K = 1, N-1
CALL DTRSM('R', 'L', 'T', 'N', NB, NB, 1D0,
& D(1,(K-1)*NB+1), LDD, B(1,(K-1)*NB+1), LDB)
CALL DSYRK('L', 'N', NB, NB, -1D0,
& B(1,(K-1)*NB+1), LDB, 1D0, D(1,K*NB+1), LDD)
CALL DPOTRF('L', NB, D(1,K*NB+1), LDD, INFO)
IF(INFO .NE. 0) THEN
INFO = INFO + K*NB
! INFO is equal to not local but global number of the row
RETURN
END IF
END DO
RETURN
END
I am using the following command to link the code to MKL and compile using the Intel visual Fortran compiler 64 bit:
python -m numpy.f2py -c --fcompiler=intelvem -L"C:\Program Files (x86)\Intel\oneAPI\compiler\2021.2.0\windows\compiler\lib\intel64_win" -lifconsol -L"C:\Program Files (x86)\Intel\oneAPI\mkl\2021.2.0\lib\intel64" -lmkl_intel_ilp64 -L"C:\Program Files (x86)\Intel\oneAPI\mkl\2021.2.0\lib\intel64" -lmkl_sequential -L"C:\Program Files (x86)\Intel\oneAPI\mkl\2021.2.0\lib\intel64" -lmkl_core -I"C:\Program Files (x86)\Intel\oneAPI\mkl\2021.2.0\include" dpbltrf.f -m SBTCF
The signature of the resulting function looks OK:
d,b,info = dpbltrf(n,nb,d,b,[ldd,ldb,overwrite_d,overwrite_b])
Wrapper for ``dpbltrf``.
Parameters
----------
n : input int
nb : input int
d : input rank-2 array('d') with bounds (ldd,*)
b : input rank-2 array('d') with bounds (ldb,*)
Other Parameters
----------------
overwrite_d : input int, optional
Default: 0
ldd : input int, optional
Default: shape(d,0)
overwrite_b : input int, optional
Default: 0
ldb : input int, optional
Default: shape(b,0)
Returns
-------
d : rank-2 array('d') with bounds (ldd,*)
b : rank-2 array('d') with bounds (ldb,*)
info : int
I can then import the created module in Python and call it. Here is a minimal example of that:
import numpy as np
import matplotlib.pyplot as plt
from SBTCF import dpbltrf
print(dpbltrf.__doc__)
def mx(x, a):
x = np.ravel(x)
npts = len(x)
distances = np.abs(x[:, None] - x[None, :])
mat = np.exp(-distances / a)
return mat
# Assemble arrays
Nx = 5
Nt = 10
t = np.linspace(0, 1, Nt)
mat_t = mx(t, 1.0)
# Stack the arrays horizontally into the shape required by dpbltrf. D is diagonal
# and B is off-diagonal
D = mat_t
for i in range(Nx-1):
D = np.hstack([D, mat_t])
B = mat_t
for i in range(Nx-2):
B = np.hstack([B, mat_t])
#plt.imshow(D)
#plt.show()
#plt.imshow(B)
#plt.show()
d, b, info = dpbltrf(Nx, Nt, D, B)
When I try to run the script I get the error:
Intel MKL ERROR: Parameter 4 was incorrect on entry to DPOTRF
I am confident that my parameters are correct and I suspect the issue is caused by wrong compilation options. Could someone tell me if there is a mistake in my compilation command that could be causing this? Or would could there be some other issue that is producing this error?
