1 /* mmm.c */
2
3 #include "../Pencil.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 --------------------------------
8 compute Y := Y + (A + sigma*B)*X
9
10 created -- 98may02, cca
11 --------------------------------
12 */
13 void
Pencil_mmm(Pencil * pencil,DenseMtx * Y,DenseMtx * X)14 Pencil_mmm (
15 Pencil *pencil,
16 DenseMtx *Y,
17 DenseMtx *X
18 ) {
19 int ncolX, ncolY, nrhs, nrow, nrowX, nrowY ;
20 /*
21 ---------------
22 check the input
23 ---------------
24 */
25 if ( pencil == NULL || Y == NULL || X == NULL ) {
26 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
27 "\n bad input\n", pencil, Y, X) ;
28 exit(-1) ;
29 }
30 if ( !(PENCIL_IS_REAL(pencil) || PENCIL_IS_COMPLEX(pencil)) ) {
31 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
32 "\n bad type %d for pencil\n", pencil, Y, X, pencil->type) ;
33 exit(-1) ;
34 }
35 if ( !(DENSEMTX_IS_REAL(Y) || DENSEMTX_IS_COMPLEX(Y)) ) {
36 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
37 "\n bad type %d for Y\n", pencil, Y, X, Y->type) ;
38 exit(-1) ;
39 }
40 if ( !(DENSEMTX_IS_REAL(X) || DENSEMTX_IS_COMPLEX(X)) ) {
41 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
42 "\n bad type %d for X\n", pencil, Y, X, X->type) ;
43 exit(-1) ;
44 }
45 if ( PENCIL_IS_REAL(pencil) && !DENSEMTX_IS_REAL(Y) ) {
46 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
47 "\n pencil is real, Y is not\n", pencil, Y, X) ;
48 exit(-1) ;
49 }
50 if ( PENCIL_IS_REAL(pencil) && !DENSEMTX_IS_REAL(X) ) {
51 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
52 "\n pencil is real, X is not\n", pencil, Y, X) ;
53 exit(-1) ;
54 }
55 if ( PENCIL_IS_COMPLEX(pencil) && !DENSEMTX_IS_COMPLEX(Y) ) {
56 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
57 "\n pencil is complex, Y is not\n", pencil, Y, X) ;
58 exit(-1) ;
59 }
60 if ( PENCIL_IS_COMPLEX(pencil) && !DENSEMTX_IS_COMPLEX(X) ) {
61 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
62 "\n pencil is complex, X is not\n", pencil, Y, X) ;
63 exit(-1) ;
64 }
65 DenseMtx_dimensions(Y, &nrowY, &ncolY) ;
66 DenseMtx_dimensions(X, &nrowX, &ncolX) ;
67 if ( nrowY != nrowX || ncolY != ncolX ) {
68 fprintf(stderr, "\n fatal error in Pencil_mmm(%p,%p,%p)"
69 "\n nrowY %d, ncolY %d, nrowX %d, ncolX %d\n",
70 pencil, Y, X, nrowY, ncolY, nrowX, ncolX) ;
71 exit(-1) ;
72 }
73 nrow = nrowY ;
74 nrhs = ncolY ;
75 if ( pencil->inpmtxA == NULL ) {
76 /*
77 -----------------
78 A is the identity
79 -----------------
80 */
81 if ( PENCIL_IS_REAL(pencil) ) {
82 double *x, *y ;
83 int irow, jrhs ;
84
85 x = DenseMtx_entries(X) ;
86 y = DenseMtx_entries(Y) ;
87 for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
88 for ( irow = 0 ; irow < nrow ; irow++ ) {
89 y[irow] += x[irow] ;
90 }
91 x += nrow ; y += nrow ;
92 }
93 } else if ( PENCIL_IS_COMPLEX(pencil) ) {
94 double *x, *y ;
95 int irow, jrhs ;
96
97 for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
98 for ( irow = 0 ; irow < nrow ; irow++ ) {
99 y[2*irow] += x[2*irow] ;
100 y[2*irow+1] += x[2*irow+1] ;
101 }
102 x += 2*nrow ; y += 2*nrow ;
103 }
104 }
105 } else {
106 double alpha[2] ;
107 /*
108 ----------------------------------------
109 A is not the identity, multiply x with A
110 ----------------------------------------
111 */
112 alpha[0] = 1.0 ; alpha[1] = 0.0 ;
113 if ( PENCIL_IS_SYMMETRIC(pencil) ) {
114 InpMtx_sym_mmm(pencil->inpmtxA, Y, alpha, X) ;
115 } else if ( PENCIL_IS_HERMITIAN(pencil) ) {
116 InpMtx_herm_mmm(pencil->inpmtxA, Y, alpha, X) ;
117 } else if ( PENCIL_IS_NONSYMMETRIC(pencil) ) {
118 InpMtx_nonsym_mmm(pencil->inpmtxA, Y, alpha, X) ;
119 }
120 }
121 if ( pencil->sigma[0] != 0.0 || pencil->sigma[1] != 0.0 ) {
122 if ( pencil->inpmtxB != NULL ) {
123 /*
124 -----------------------------------------
125 B is not the identity, add sigma*B*x to y
126 -----------------------------------------
127 */
128 if ( PENCIL_IS_SYMMETRIC(pencil) ) {
129 InpMtx_sym_mmm(pencil->inpmtxB, Y, pencil->sigma, X) ;
130 } else if ( PENCIL_IS_HERMITIAN(pencil) ) {
131 InpMtx_herm_mmm(pencil->inpmtxB, Y, pencil->sigma, X) ;
132 } else if ( PENCIL_IS_NONSYMMETRIC(pencil) ) {
133 InpMtx_nonsym_mmm(pencil->inpmtxB, Y, pencil->sigma, X) ;
134 }
135 } else {
136 /*
137 -----------------------------------
138 B is the identity, add sigma*x to y
139 -----------------------------------
140 */
141 if ( PENCIL_IS_REAL(pencil) ) {
142 double sigmareal ;
143 double *x, *y ;
144 int irow, jrhs ;
145
146 x = DenseMtx_entries(X) ; y = DenseMtx_entries(Y) ;
147 sigmareal = pencil->sigma[0] ;
148 for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
149 for ( irow = 0 ; irow < nrow ; irow++ ) {
150 y[irow] += sigmareal*x[irow] ;
151 }
152 x += nrow ; y += nrow ;
153 }
154 } else if ( PENCIL_IS_COMPLEX(pencil) ) {
155 double sigmaimag, sigmareal, ximag, xreal ;
156 double *x, *y ;
157 int irow, jrhs ;
158
159 x = DenseMtx_entries(X) ; y = DenseMtx_entries(Y) ;
160 sigmareal = pencil->sigma[0] ; sigmaimag = pencil->sigma[1] ;
161 for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
162 for ( irow = 0 ; irow < nrow ; irow++ ) {
163 xreal = x[2*irow] ; ximag = x[2*irow+1] ;
164 y[2*irow] += sigmareal*xreal - sigmaimag*ximag ;
165 y[2*irow+1] += sigmareal*ximag + sigmaimag*xreal ;
166 }
167 x += 2*nrow ; y += 2*nrow ;
168 }
169 }
170 }
171 }
172 return ; }
173
174 /*--------------------------------------------------------------------*/
175