160
161
162
163
164
165 IMPLICIT NONE
166
167
168 CHARACTER SIDE, TRANS
169 INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
170
171
172 COMPLEX Q( LDQ, * ), C( LDC, * ), WORK( * )
173
174
175
176
177
178 COMPLEX ONE
179 parameter( one = ( 1.0e+0, 0.0e+0 ) )
180
181
182 LOGICAL LEFT, LQUERY, NOTRAN
183 INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
184
185
186 LOGICAL LSAME
188
189
191
192
193 INTRINSIC cmplx, max, min
194
195
196
197
198
199 info = 0
200 left =
lsame( side,
'L' )
201 notran =
lsame( trans,
'N' )
202 lquery = ( lwork.EQ.-1 )
203
204
205
206
207 IF( left ) THEN
208 nq = m
209 ELSE
210 nq = n
211 END IF
212 nw = nq
213 IF( n1.EQ.0 .OR. n2.EQ.0 ) nw = 1
214 IF( .NOT.left .AND. .NOT.
lsame( side,
'R' ) )
THEN
215 info = -1
216 ELSE IF( .NOT.
lsame( trans,
'N' ) .AND.
217 $ .NOT.
lsame( trans,
'C' ) )
218 $ THEN
219 info = -2
220 ELSE IF( m.LT.0 ) THEN
221 info = -3
222 ELSE IF( n.LT.0 ) THEN
223 info = -4
224 ELSE IF( n1.LT.0 .OR. n1+n2.NE.nq ) THEN
225 info = -5
226 ELSE IF( n2.LT.0 ) THEN
227 info = -6
228 ELSE IF( ldq.LT.max( 1, nq ) ) THEN
229 info = -8
230 ELSE IF( ldc.LT.max( 1, m ) ) THEN
231 info = -10
232 ELSE IF( lwork.LT.nw .AND. .NOT.lquery ) THEN
233 info = -12
234 END IF
235
236 IF( info.EQ.0 ) THEN
237 lwkopt = m*n
238 work( 1 ) = cmplx( lwkopt )
239 END IF
240
241 IF( info.NE.0 ) THEN
242 CALL xerbla(
'CUNM22', -info )
243 RETURN
244 ELSE IF( lquery ) THEN
245 RETURN
246 END IF
247
248
249
250 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
251 work( 1 ) = 1
252 RETURN
253 END IF
254
255
256
257 IF( n1.EQ.0 ) THEN
258 CALL ctrmm( side,
'Upper', trans,
'Non-Unit', m, n, one,
259 $ q, ldq, c, ldc )
260 work( 1 ) = one
261 RETURN
262 ELSE IF( n2.EQ.0 ) THEN
263 CALL ctrmm( side,
'Lower', trans,
'Non-Unit', m, n, one,
264 $ q, ldq, c, ldc )
265 work( 1 ) = one
266 RETURN
267 END IF
268
269
270
271 nb = max( 1, min( lwork, lwkopt ) / nq )
272
273 IF( left ) THEN
274 IF( notran ) THEN
275 DO i = 1, n, nb
276 len = min( nb, n-i+1 )
277 ldwork = m
278
279
280
281 CALL clacpy(
'All', n1, len, c( n2+1, i ), ldc, work,
282 $ ldwork )
283 CALL ctrmm(
'Left',
'Lower',
'No Transpose',
284 $ 'Non-Unit',
285 $ n1, len, one, q( 1, n2+1 ), ldq, work,
286 $ ldwork )
287
288
289
290 CALL cgemm(
'No Transpose',
'No Transpose', n1, len,
291 $ n2,
292 $ one, q, ldq, c( 1, i ), ldc, one, work,
293 $ ldwork )
294
295
296
297 CALL clacpy(
'All', n2, len, c( 1, i ), ldc,
298 $ work( n1+1 ), ldwork )
299 CALL ctrmm(
'Left',
'Upper',
'No Transpose',
300 $ 'Non-Unit',
301 $ n2, len, one, q( n1+1, 1 ), ldq,
302 $ work( n1+1 ), ldwork )
303
304
305
306 CALL cgemm(
'No Transpose',
'No Transpose', n2, len,
307 $ n1,
308 $ one, q( n1+1, n2+1 ), ldq, c( n2+1, i ), ldc,
309 $ one, work( n1+1 ), ldwork )
310
311
312
313 CALL clacpy(
'All', m, len, work, ldwork, c( 1, i ),
314 $ ldc )
315 END DO
316 ELSE
317 DO i = 1, n, nb
318 len = min( nb, n-i+1 )
319 ldwork = m
320
321
322
323 CALL clacpy(
'All', n2, len, c( n1+1, i ), ldc, work,
324 $ ldwork )
325 CALL ctrmm(
'Left',
'Upper',
'Conjugate',
'Non-Unit',
326 $ n2, len, one, q( n1+1, 1 ), ldq, work,
327 $ ldwork )
328
329
330
331 CALL cgemm(
'Conjugate',
'No Transpose', n2, len, n1,
332 $ one, q, ldq, c( 1, i ), ldc, one, work,
333 $ ldwork )
334
335
336
337 CALL clacpy(
'All', n1, len, c( 1, i ), ldc,
338 $ work( n2+1 ), ldwork )
339 CALL ctrmm(
'Left',
'Lower',
'Conjugate',
'Non-Unit',
340 $ n1, len, one, q( 1, n2+1 ), ldq,
341 $ work( n2+1 ), ldwork )
342
343
344
345 CALL cgemm(
'Conjugate',
'No Transpose', n1, len, n2,
346 $ one, q( n1+1, n2+1 ), ldq, c( n1+1, i ), ldc,
347 $ one, work( n2+1 ), ldwork )
348
349
350
351 CALL clacpy(
'All', m, len, work, ldwork, c( 1, i ),
352 $ ldc )
353 END DO
354 END IF
355 ELSE
356 IF( notran ) THEN
357 DO i = 1, m, nb
358 len = min( nb, m-i+1 )
359 ldwork = len
360
361
362
363 CALL clacpy(
'All', len, n2, c( i, n1+1 ), ldc, work,
364 $ ldwork )
365 CALL ctrmm(
'Right',
'Upper',
'No Transpose',
366 $ 'Non-Unit',
367 $ len, n2, one, q( n1+1, 1 ), ldq, work,
368 $ ldwork )
369
370
371
372 CALL cgemm(
'No Transpose',
'No Transpose', len, n2,
373 $ n1,
374 $ one, c( i, 1 ), ldc, q, ldq, one, work,
375 $ ldwork )
376
377
378
379 CALL clacpy(
'All', len, n1, c( i, 1 ), ldc,
380 $ work( 1 + n2*ldwork ), ldwork )
381 CALL ctrmm(
'Right',
'Lower',
'No Transpose',
382 $ 'Non-Unit',
383 $ len, n1, one, q( 1, n2+1 ), ldq,
384 $ work( 1 + n2*ldwork ), ldwork )
385
386
387
388 CALL cgemm(
'No Transpose',
'No Transpose', len, n1,
389 $ n2,
390 $ one, c( i, n1+1 ), ldc, q( n1+1, n2+1 ), ldq,
391 $ one, work( 1 + n2*ldwork ), ldwork )
392
393
394
395 CALL clacpy(
'All', len, n, work, ldwork, c( i, 1 ),
396 $ ldc )
397 END DO
398 ELSE
399 DO i = 1, m, nb
400 len = min( nb, m-i+1 )
401 ldwork = len
402
403
404
405 CALL clacpy(
'All', len, n1, c( i, n2+1 ), ldc, work,
406 $ ldwork )
407 CALL ctrmm(
'Right',
'Lower',
'Conjugate',
'Non-Unit',
408 $ len, n1, one, q( 1, n2+1 ), ldq, work,
409 $ ldwork )
410
411
412
413 CALL cgemm(
'No Transpose',
'Conjugate', len, n1, n2,
414 $ one, c( i, 1 ), ldc, q, ldq, one, work,
415 $ ldwork )
416
417
418
419 CALL clacpy(
'All', len, n2, c( i, 1 ), ldc,
420 $ work( 1 + n1*ldwork ), ldwork )
421 CALL ctrmm(
'Right',
'Upper',
'Conjugate',
'Non-Unit',
422 $ len, n2, one, q( n1+1, 1 ), ldq,
423 $ work( 1 + n1*ldwork ), ldwork )
424
425
426
427 CALL cgemm(
'No Transpose',
'Conjugate', len, n2, n1,
428 $ one, c( i, n2+1 ), ldc, q( n1+1, n2+1 ), ldq,
429 $ one, work( 1 + n1*ldwork ), ldwork )
430
431
432
433 CALL clacpy(
'All', len, n, work, ldwork, c( i, 1 ),
434 $ ldc )
435 END DO
436 END IF
437 END IF
438
439 work( 1 ) = cmplx( lwkopt )
440 RETURN
441
442
443
subroutine xerbla(srname, info)
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
logical function lsame(ca, cb)
LSAME
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM