home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / comp / lang / c / 17010 < prev    next >
Encoding:
Internet Message Format  |  1992-11-22  |  3.2 KB
  1. From: gajdos@hpavla.lf.hp.com (Larry Gajdos)
  2. Date: Thu, 19 Nov 1992 15:13:40 GMT
  3. Subject: Re: Reasons for using C vs. Fortran or vice/versa
  4. Message-ID: <9130057@hpavla.lf.hp.com>
  5. Organization: Hewlett-Packard Little Falls Site
  6. Path: sparky!uunet!cs.utexas.edu!sdd.hp.com!apollo.hp.com!cupnews0.cup.hp.com!nsa.hp.com!hpscit.sc.hp.com!scd.hp.com!hpscdm!hplextra!hpcss01!hpwala!hpavla!gajdos
  7. Newsgroups: comp.lang.c
  8. References: <1992Nov12.135901.15191@ccd.harris.com>
  9. Lines: 87
  10.  
  11. In comp.lang.c, rons@hardy.u.washington.edu (Ronald Schoenberg) writes:
  12.  
  13. > In article <gay.721614205@sfu.ca> gay@selkirk.sfu.ca (Ian D. Gay) writes:
  14. > >rons@hardy.u.washington.edu (Ronald Schoenberg) writes:
  15. > >
  16. > >>  Fortran arrays are column echelon for a
  17. > >>reason, and C arrays are row echelon slowing down matrix operations
  18. > >>relative to Fortran.     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  19. > >
  20. > >Maybe I'm missing something, but I don't see why this should be so. 
  21. > >Could you explain more fully?
  22. > >
  23. > >
  24. > >>Ron Schoenberg
  25.  
  26. > The most important matrix operation in numerical work is the gaxpy
  27. > ("generalized a times x plus y").  Numerical routines are written to
  28. > be rich in gaxpies (see Matrix Computations by Golub and Van Loan for
  29. > details). In the source code below I have written two gaxpies, one
  30. > for data stored column echelon and a second for the same data stored
  31. > row echelon.  The only difference between the two is in the
  32. > calculation of the index of the element of the two dimensional array
  33. > used in the multiplication.  You will note that the index in the
  34. > column echelon storage is a simple increment, whereas in the row
  35. > echelon storage two operations in addition to the assignment are 
  36. > required.  I know very little about assembly code, but I believe that an
  37. > assembly version of the column echelon gaxpy is especially efficient,
  38. > while the assembly code for row echelon would involve a number of extra
  39. > operations.
  40.  
  41. A simple rearrangment of the row echelon calculation is shown below. This
  42. seems to show that there is no difference! Or am I missing something?
  43.  
  44. > #include <stdio.h>
  45.  
  46. > void main()
  47. > {
  48. >     float a_col[] = { 1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16 };
  49. >     float a_row[] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 };
  50. >     float x[] = { 1,2,3,4 };
  51. >     float z[4];
  52.  
  53. >     int m = 4, n = 4;
  54. >     int i, j, k;
  55.  
  56. > /*
  57. > ** column echelon gaxpy
  58. > */
  59. >     for (i = 0; i < m; i++) {
  60. >         z[i] = 0;
  61. >         k = i;
  62. >         for (j = 0; j < n; j++) {
  63. >             z[i] += a_col[k] * x[j];
  64. >             k += n;
  65. >         }
  66. >     }
  67.  
  68. > /*
  69. > ** row echelon gaxpy
  70. > */
  71. >     for (i = 0; i < m; i++) {  
  72. >         z[i] = 0;
  73. >         for (j = 0; j < n; j++) {
  74. >             k = i * n + j;
  75. >             z[i] += a_row[k] * x[j];
  76. >         }
  77. >     }
  78.  
  79.  /***** Redo the above row echelon gaxpy: *****/
  80.  
  81.       k = 0;
  82.       for (i = 0; i < m; i++) {  
  83.           z[i] = 0;
  84.           /***** k will be n*i on entry to loop *****/
  85.           for (j = 0; j < n; j++) {
  86.               /***** k will be n*i + j on each pass thru loop *****/
  87.               z[i] += a_row[k] * x[j];
  88.               k++;
  89.           }
  90.       }
  91.  
  92. > }
  93.  
  94.  
  95. If anything, the code will be faster for the modfied row echelon calculation.
  96.  
  97. Larry Gajdos
  98.