A gentle critque..

[Paulo Herrera]

Walter Bright wrote:
> Paulo Herrera wrote:
>> Walter Bright wrote:
>>> Don Clugston wrote:
>>>
>>>>> About native libraries
>>>>> ----------------------
>>>>> I think the way to go is to create QUALITY native D libraries. I 
>>>>> think we usually overestimate the cost of porting/developing native 
>>>>> libraries.
>>>>
>>>> I agree. When porting a library
>>>> (a) there isn't any algorithm development (not much 'thinking' time 
>>>> required);
>>>> (b) it's very easy to test (you can run test cases against the 
>>>> original version);
>>>> (c) the D version is frequently greatly superior to the original.
>>>
>>> Hmm. Is there a 'canonical' Fortran numerics library with a solid 
>>> test suite that could be translated to D as a showcase project?
>>
>> Do you have some idea in mind?
> 
> I am not familiar with the various Fortran libraries out there, so I 
> have no specific one in mind.
> 
> 
>> If we can make a port of that library to D and show that it performs 
>> close to  SPARSKIT, that would be a good demonstration of D 
>> capabilities for numerical computing.
> 
> Want to get started on it? <g>
OK, I take the challenge. However, I'm pretty new to D so I'd like to 
ask some questions before starting.

As someone else posted, most of numerical libraries in Fortran and C/C++ 
are based on BLAS and LAPACK. So, a first logic step is to evaluate if 
we should port those libraries to D. This weekend I took a look to the 
specification of those libraries (http://www.netlib.org/blas/). After 
going through the document I'm not sure what is the best way to write a 
similar library in a "D way". I will try to explain what I mean by "D 
way" ....
I'm sure many people are familiar with those libraries, but I include a 
brief explanation below to make my question clear.

Those libraries define functions to compute basic matrix/vector 
operations, such as:
- x = y, where x and y are vectors
- R = aAx, where A and R are matrices, a is a scalar, and x is a vector
- r = aAx + by; where a and b are scalars, A is a matrix, and r, x and y 
  are vectors.

Since they were originally developed as a specification for old 
Fortran77 and C, the declaration of the routines that implement those 
operations look like:
- x = y    => void copy(T *x, T *y)
- R = aAx  => void mult0(T a, T *A, T *x, T *r), etc.

Those declarations are not elegant and they are not easy to use. My 
first question is: Should we write a D library in that way to get best 
performance? I hope the answer is no, because that would be an advantage 
of D over other languages.
I believe we should overwrite operators to make the notation as natural 
as possible. This is something that the designers of BLAS also faced 
about Fortran95 that includes array operations: "Some of the functions 
... can be replaced by simple array expressions and assignments in 
Fortran95, without loss of convenience or performance(assuming a 
reasonable degree of optimization by the compiler)...." (in "Basic 
Linear Algebra Subprograms Technical (BLAST) Forum Standard", 2001, p. 26)

My questions are:
1) Do you think it's possible to write vector and matrix classes with 
overloaded operators that perform as well as the primitive BLAS 
operations? What about temporary objects?
2) Do you think some of those operations can be accelerated by 
implementing them as part of the language (as described in future work 
in D) If yes, is there any time frame for those changes?
3) Since the type T on those expressions can change, what is the best 
way to implement those function without loosing performance?  templates? 
   Since in general people will use only one data type in their 
programs, could we use typedef or alias to get better performance?
4) Is there a real interest for this kind of libraries?

Thanks,
Paulo.