home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / eigen / blas.c next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.5 KB  |  158 lines
  1. /* Assorted matrix functions.
  2.  */
  3.  
  4.  
  5. /* Multiply r (rows) by c (columns) matrix A on the left
  6.  * by column vector V of dimension c on the right
  7.  * to produce a (column) vector Y output of dimension r.
  8.  */
  9. mvmpy( r, c, A, V, Y )
  10. int r, c;
  11. double *A, *V, *Y;
  12. {
  13. register double s;
  14. double *pA, *pV, *pY;
  15. int i, j;
  16.  
  17. pA = A;
  18. pY = Y;
  19. for( i=0; i<r; i++ )
  20.     {
  21.     pV = V;
  22.     s = 0.0;
  23.     for( j=0; j<c; j++ )
  24.         {
  25.         s += *pA++ * *pV++;
  26.         }
  27.     *pY++ = s;
  28.     }
  29. }
  30.  
  31.  
  32. /* Multiply an r (rows) by c (columns) matrix A on the left
  33.  * by a c (rows) by r (columns) matrix B on the right
  34.  * to produce an r by r matrix Y.
  35.  */
  36. mmmpy( r, c, A, B, Y )
  37. int r, c;
  38. double *A, *B, *Y;
  39. {
  40. register double s;
  41. double *pA, *pB, *pY, *pt;
  42. int i, j, k;
  43.  
  44. pY = Y;
  45. pB = B;
  46. for( i=0; i<r; i++ )
  47.     {
  48.     pA = A;
  49.     for( j=0; j<r; j++ )
  50.         {
  51.         pt = pB;
  52.         s = 0.0;
  53.         for( k=0; k<c; k++ )
  54.             {
  55.             s += *pA++ * *pt;
  56.             pt += r; /* increment to next row underneath */
  57.             }
  58.         *pY++ = s;
  59.         }
  60.     pB += 1;
  61.     }
  62. }
  63.  
  64.  
  65. /* Transpose the n by n square matrix A and put the result in T.
  66.  * T may occupy the same storage as A.
  67.  */
  68. mtransp( n, A, T )
  69. int n;
  70. double *A, *T;
  71. {
  72. int i, j, np1;
  73. double *pAc, *pAr, *pTc, *pTr, *pA0, *pT0;
  74. double x, y;
  75.  
  76. np1 = n+1;
  77. pA0 = A;
  78. pT0 = T;
  79. for( i=0; i<n-1; i++ ) /* row index */
  80.     {
  81.     pAc = pA0; /* next diagonal element of input */
  82.     pAr = pAc + n; /* next row down underneath the diagonal element */
  83.     pTc = pT0; /* next diagonal element of the output */
  84.     pTr = pTc + n; /* next row underneath */
  85.     *pTc++ = *pAc++; /* copy the diagonal element */
  86.     for( j=i+1; j<n; j++ ) /* column index */
  87.         {
  88.         x = *pAr;
  89.         *pTr = *pAc++;
  90.         *pTc++ = x;
  91.         pAr += n;
  92.         pTr += n;
  93.         }
  94.     pA0 += np1; /* &A[n*i+i] for next i */
  95.     pT0 += np1; /* &T[n*i+i] for next i */
  96.     }
  97. *pT0 = *pA0; /* copy the diagonal element */
  98. }
  99.  
  100.  
  101. /* Return maximum off-diagonal element of n by n square matrix A
  102.  */
  103. double maxoffd( n, A )
  104. int n;
  105. double *A;
  106. {
  107. double e, x;
  108. int i, j, nm1;
  109. double *pA;
  110.  
  111. nm1 = n-1;
  112. e = 0.0;
  113. pA = A;
  114. for( i=0; i<nm1; i++ )
  115.     {
  116.     ++pA; /* skip over the diagonal element */
  117.     for( j=0; j<n; j++ )
  118.         {
  119.         x = *pA++;
  120.         if( x < 0 )
  121.             x = -x;
  122.         if( x > e )
  123.             e = x;
  124.         }
  125.     }
  126. return( e );
  127. }
  128.  
  129.  
  130.  
  131.  
  132. /* Unpack symmetric matrix T stored in lower triangular form
  133.  * into a symmetric n by n square matrix S.
  134.  */
  135. tritosquare( n, T, S )
  136. int n;
  137. double T[], S[];
  138. {
  139. double *pT;
  140. int i, j, ni, nj;
  141.  
  142. /* offset to (i,j) element is (j*j+j)/2 + i */
  143. pT = T;
  144. ni = 0;
  145. for( i=0; i<n; i++ )
  146.     {
  147.     nj = 0;
  148.     for( j=0; j<i; j++ )
  149.         {
  150.         S[ni+j] = *pT;
  151.         S[nj+i] = *pT++;
  152.         nj += n;
  153.         }
  154.     S[ni+i] = *pT++;
  155.     ni += n;
  156.     }
  157. }
  158.