home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / comp / lang / c / 16445 < prev    next >
Encoding:
Text File  |  1992-11-13  |  3.4 KB  |  95 lines
  1. Path: sparky!uunet!europa.asd.contel.com!darwin.sura.net!wupost!uwm.edu!ogicse!news.u.washington.edu!hardy.u.washington.edu!rons
  2. From: rons@hardy.u.washington.edu (Ronald Schoenberg)
  3. Newsgroups: comp.lang.c
  4. Subject: Re: Reasons for using C vs. Fortran or vice/versa
  5. Message-ID: <1992Nov13.060701.6917@u.washington.edu>
  6. Date: 13 Nov 92 06:07:01 GMT
  7. Article-I.D.: u.1992Nov13.060701.6917
  8. References: <1992Nov12.135901.15191@ccd.harris.com> <1992Nov12.205823.13951@u.washington.edu> <gay.721614205@sfu.ca>
  9. Sender: news@u.washington.edu (USENET News System)
  10. Distribution: usa
  11. Organization: University of Washington, Seattle
  12. Lines: 81
  13.  
  14. In article <gay.721614205@sfu.ca> gay@selkirk.sfu.ca (Ian D. Gay) writes:
  15. >rons@hardy.u.washington.edu (Ronald Schoenberg) writes:
  16. >
  17. >>  Fortran arrays are column echelon for a
  18. >>reason, and C arrays are row echelon slowing down matrix operations
  19. >>relative to Fortran.     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  20. >
  21. >Maybe I'm missing something, but I don't see why this should be so. 
  22. >Could you explain more fully?
  23. >
  24. >
  25. >>Ron Schoenberg
  26.  
  27. The most important matrix operation in numerical work is the gaxpy
  28. ("generalized a times x plus y").  Numerical routines are written to
  29. be rich in gaxpies (see Matrix Computations by Golub and Van Loan for
  30. details). In the source code below I have written two gaxpies, one
  31. for data stored column echelon and a second for the same data stored
  32. row echelon.  The only difference between the two is in the
  33. calculation of the index of the element of the two dimensional array
  34. used in the multiplication.  You will note that the index in the
  35. column echelon storage is a simple increment, whereas in the row
  36. echelon storage two operations in addition to the assignment are 
  37. required.  I know very little about assembly code, but I believe that an
  38. assembly version of the column echelon gaxpy is especially efficient,
  39. while the assembly code for row echelon would involve a number of extra
  40. operations.
  41.  
  42. As I'm sure you well know, if you were to actually use the two
  43. dimensional arrays in C (I avoid it myself), the generated 
  44. assembly code would contain the extra operations because data in
  45. two-dimensional C arrays are row echelon.  C programmers have
  46. dispensed with one additional operation by using the zero offset.
  47. But for really efficient numerical work that is not enough.
  48.  
  49. Actually I program in C++ using the M++ matrix class library.  
  50. M++ is based on the BLAS (Basic Linear Algebra Subroutines) which are
  51. most efficient when data are stored column echelon, and therefore
  52. M++ stores its data column echelon.  But I don't have
  53. to think very much about how the data are stored because that is all
  54. taken care of for me behind the scenes (viva la C++).
  55.  
  56. I apologize for any faux pas in the code below - I am a converted
  57. Fortran number-cruncher.
  58.  
  59.  
  60. #include <stdio.h>
  61.  
  62. void main()
  63. {
  64.     float a_col[] = { 1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16 };
  65.     float a_row[] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 };
  66.     float x[] = { 1,2,3,4 };
  67.     float z[4];
  68.  
  69.     int m = 4, n = 4;
  70.     int i, j, k;
  71.  
  72. /*
  73. ** column echelon gaxpy
  74. */
  75.     for (i = 0; i < m; i++) {
  76.         z[i] = 0;
  77.         k = i;
  78.         for (j = 0; j < n; j++) {
  79.             z[i] += a_col[k] * x[j];
  80.             k += n;
  81.         }
  82.     }
  83.  
  84. /*
  85. ** row echelon gaxpy
  86. */
  87.     for (i = 0; i < m; i++) {  
  88.         z[i] = 0;
  89.         for (j = 0; j < n; j++) {
  90.             k = i * n + j;
  91.             z[i] += a_row[k] * x[j];
  92.         }
  93.     }
  94. }
  95.