Skip to content

Is unused triangular in DPOTRF guranteed to be 0? #1408

Description

@dasergatskov

For DPOTRF('U',...), does the documented statement that the strictly lower triangular part of A “is not referenced” guarantee that it is not modified, or may an implementation use it as scratch storage even though only the upper factor is specified on output?
Octave's code seems to assume so and that worked fine with ref. Lapack and with OpenBlas, but breaks with Apple's both Veclib and ACCELERATE_NEW_LAPACK.
Here is the test F77 program:

chol33.f.gz

 gfortran -framework Accelerate chol33.f  

and

 % ./a.out 
 DPOTRF INFO =    0

 Columns  1 through  8:

    1.4142    0.7071    0.7071    0.7071    0.7071    0.7071    0.7071    0.7071
   -1.0000    1.2247    0.4082    0.4082    0.4082    0.4082    0.4082    0.4082
   -0.5918   -1.0000    1.1547    0.2887    0.2887    0.2887    0.2887    0.2887
   -0.4209   -0.7113   -1.0000    1.1180    0.2236    0.2236    0.2236    0.2236
   -0.3268   -0.5523   -0.7764   -1.0000    1.0954    0.1826    0.1826    0.1826
   -0.2671   -0.4514   -0.6346   -0.8174   -1.0000    1.0801    0.1543    0.1543
   -0.2259   -0.3818   -0.5367   -0.6913   -0.8457   -1.0000    1.0690    0.1336
   -0.1957   -0.3308   -0.4650   -0.5989   -0.7327   -0.8664   -1.0000    1.0607
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
...
 Nonzeros in strict lower triangle:   56 of  528
 R is NOT upper triangular

Linking to reference lapack produces "expected" results:

% gfortran -L/opt/homebrew/Cellar/lapack/3.12.1_1/lib -llapack chol33.f
% ./a.out                                                              
 DPOTRF INFO =    0

 Columns  1 through  8:

    1.4142    0.7071    0.7071    0.7071    0.7071    0.7071    0.7071    0.7071
    0.0000    1.2247    0.4082    0.4082    0.4082    0.4082    0.4082    0.4082
    0.0000    0.0000    1.1547    0.2887    0.2887    0.2887    0.2887    0.2887
...
 Nonzeros in strict lower triangle:    0 of  528
 R is upper triangular

So, is it an Apple's bug or our expectations are wrong?

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions