MAGMA  1.2.0
MatrixAlgebraonGPUandMulticoreArchitectures
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
zgeqrf_mc.cpp
Go to the documentation of this file.
1 /*
2  -- MAGMA (version 1.2.0) --
3  Univ. of Tennessee, Knoxville
4  Univ. of California, Berkeley
5  Univ. of Colorado, Denver
6  May 2012
7 
8  @precisions normal z -> s d c
9 
10 */
11 #include "common_magma.h"
12 
13 #define A(m,n) (a+(n)*(*lda)+(m))
14 #define T(m) (work+(m)*(nb))
15 #define W(k,n) &(local_work[(mt)*(n-1)+(k)])
16 
17 void getro (char *trans, const magma_int_t m, const magma_int_t n,
18  const cuDoubleComplex *A, const magma_int_t LDA,
19  cuDoubleComplex *B, const magma_int_t LDB)
20 {
21  const cuDoubleComplex *Atmp;
22 
23  magma_int_t i, j;
24 
25  for (i=0; i<m; i++) {
26  Atmp = A + i;
27  for (j=0; j<n; j++) {
28 
29  if (trans[0] == 'C') {
30 
31  double *ap,*tp;
32  cuDoubleComplex tmp;
33  ap = (double *)&(*Atmp);
34  tp = (double *)&tmp;
35  tp[0] = ap[0];
36  tp[1] = -ap[1];
37 
38  *B = tmp;
39 
40  } else {
41 
42  *B = *Atmp;
43 
44  }
45 
46  B += 1;
47  Atmp += LDA;
48  }
49  B += (LDB - n);
50  }
51 }
52 
53 // task execution code
54 void SCHED_zlarfb(Quark* quark)
55 {
56  magma_int_t M;
57  magma_int_t N;
58  magma_int_t MM;
59  magma_int_t NN;
60  magma_int_t IB;
61  magma_int_t LDV;
62  magma_int_t LDC;
63  magma_int_t LDT;
64  magma_int_t LDW;
65 
66  cuDoubleComplex *V;
67  cuDoubleComplex *C;
68  cuDoubleComplex *T;
69  cuDoubleComplex **W;
70 
71  quark_unpack_args_13(quark, M, N, MM, NN, IB, V, LDV, C, LDC, T, LDT, W, LDW);
72 
73  if (M < 0) {
74  printf("SCHED_zlarfb: illegal value of M\n");
75  }
76  if (N < 0) {
77  printf("SCHED_zlarfb: illegal value of N\n");
78  }
79  if (IB < 0) {
80  printf("SCHED_zlarfb: illegal value of IB\n");
81  }
82 
83  *W = (cuDoubleComplex*) malloc(LDW*MM*sizeof(cuDoubleComplex));
84 
85  getro(MagmaConjTransStr, MM, NN, C, LDC, *W, LDW);
86 
87  cuDoubleComplex c_one = MAGMA_Z_ONE;
88 
89  blasf77_ztrmm("r","l","n","u",
90  &NN, &MM, &c_one, V, &LDV, *W, &LDW);
91 
92  magma_int_t K=M-MM;
93 
94  blasf77_zgemm(MagmaConjTransStr, "n", &NN, &MM, &K,
95  &c_one, &C[MM], &LDC, &V[MM], &LDV, &c_one, *W, &LDW);
96 
97  blasf77_ztrmm("r", "u", "n", "n",
98  &NN, &MM, &c_one, T, &LDT, *W, &LDW);
99 
100 }
101 
102 // task execution code
103 void SCHED_zgeqrt(Quark* quark)
104 {
105  magma_int_t M;
106  magma_int_t N;
107  magma_int_t IB;
108  cuDoubleComplex *A;
109  magma_int_t LDA;
110  cuDoubleComplex *T;
111  magma_int_t LDT;
112  cuDoubleComplex *TAU;
113  cuDoubleComplex *WORK;
114 
115  magma_int_t iinfo;
116  magma_int_t lwork=-1;
117 
118  quark_unpack_args_9(quark, M, N, IB, A, LDA, T, LDT, TAU, WORK);
119 
120  if (M < 0) {
121  printf("SCHED_zgeqrt: illegal value of M\n");
122  }
123 
124  if (N < 0) {
125  printf("SCHED_zgeqrt: illegal value of N\n");
126  }
127 
128  if ((IB < 0) || ( (IB == 0) && ((M > 0) && (N > 0)) )) {
129  printf("SCHED_zgeqrt: illegal value of IB\n");
130  }
131 
132  if ((LDA < max(1,M)) && (M > 0)) {
133  printf("SCHED_zgeqrt: illegal value of LDA\n");
134  }
135 
136  if ((LDT < max(1,IB)) && (IB > 0)) {
137  printf("SCHED_zgeqrt: illegal value of LDT\n");
138  }
139 
140  lapackf77_zgeqrf(&M, &N, A, &LDA, TAU, WORK, &lwork, &iinfo);
141  lwork=(magma_int_t)MAGMA_Z_REAL(WORK[0]);
142  lapackf77_zgeqrf(&M, &N, A, &LDA, TAU, WORK, &lwork, &iinfo);
143 
144  lapackf77_zlarft("F", "C", &M, &N, A, &LDA, TAU, T, &LDT);
145 
146 }
147 
148 // task execution code
149 void SCHED_ztrmm(Quark *quark)
150 {
151  magma_int_t m;
152  magma_int_t n;
153  cuDoubleComplex alpha;
154  cuDoubleComplex *a;
155  magma_int_t lda;
156  cuDoubleComplex **b;
157  magma_int_t ldb;
158  cuDoubleComplex beta;
159  cuDoubleComplex *c;
160  magma_int_t ldc;
161  cuDoubleComplex *work;
162 
163  magma_int_t j;
164 
165  magma_int_t one = 1;
166 
167  quark_unpack_args_11(quark, m, n, alpha, a, lda, b, ldb, beta, c, ldc, work);
168 
169  if (m < 0) {
170  printf("SCHED_ztrmm: illegal value of m\n");
171  }
172 
173  if (n < 0) {
174  printf("SCHED_ztrmm: illegal value of n\n");
175  }
176 
177  getro(MagmaConjTransStr, n, m, *b, ldb, work, m);
178 
179  blasf77_ztrmm("l", "l", "n", "u",
180  &m, &n, &alpha, a, &lda, work, &m);
181 
182  for (j = 0; j < n; j++)
183  {
184  blasf77_zaxpy(&m, &beta, &(work[j*m]), &one, &(c[j*ldc]), &one);
185  }
186 
187 }
188 
189 // task execution code
190 static void SCHED_zgemm(Quark *quark)
191 {
192  magma_int_t m;
193  magma_int_t n;
194  magma_int_t k;
195  cuDoubleComplex alpha;
196  cuDoubleComplex *a;
197  magma_int_t lda;
198  cuDoubleComplex **b;
199  magma_int_t ldb;
200  cuDoubleComplex beta;
201  cuDoubleComplex *c;
202  magma_int_t ldc;
203 
204  cuDoubleComplex *fake;
205 
206  magma_int_t dkdk;
207 
208  quark_unpack_args_13(quark, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc, fake, dkdk);
209 
210  blasf77_zgemm("n", MagmaConjTransStr,
211  &m, &n, &k, &alpha, a, &lda, *b, &ldb, &beta, c, &ldc);
212 
213 }
214 
215 // wrapper around QUARK_Insert_Task .. makes code more readable at task insertion time
216 void QUARK_Insert_Task_zgemm(Quark *quark, Quark_Task_Flags *task_flags,
217  magma_int_t m,
218  magma_int_t n,
219  magma_int_t k,
220  cuDoubleComplex alpha,
221  cuDoubleComplex *a,
222  magma_int_t lda,
223  cuDoubleComplex **b,
224  magma_int_t ldb,
225  cuDoubleComplex beta,
226  cuDoubleComplex *c,
227  magma_int_t ldc,
228  cuDoubleComplex *fake,
229  char *dag_label,
230  magma_int_t priority,
231  magma_int_t dkdk)
232 {
233 
234  QUARK_Insert_Task(quark, SCHED_zgemm, task_flags,
235  sizeof(magma_int_t), &m, VALUE,
236  sizeof(magma_int_t), &n, VALUE,
237  sizeof(magma_int_t), &k, VALUE,
238  sizeof(cuDoubleComplex), &alpha, VALUE,
239  sizeof(cuDoubleComplex)*ldb*ldb, a, INPUT,
240  sizeof(magma_int_t), &lda, VALUE,
241  sizeof(cuDoubleComplex*), b, INPUT,
242  sizeof(magma_int_t), &ldb, VALUE,
243  sizeof(cuDoubleComplex), &beta, VALUE,
244  sizeof(cuDoubleComplex)*ldb*ldb, c, INOUT | LOCALITY,
245  sizeof(magma_int_t), &ldc, VALUE,
246  sizeof(cuDoubleComplex)*ldb*ldb, fake, OUTPUT | GATHERV,
247  sizeof(magma_int_t), &priority, VALUE | TASK_PRIORITY,
248  sizeof(magma_int_t),&dkdk,VALUE,
249  strlen(dag_label)+1, dag_label, VALUE | TASKLABEL,
250 6, "purple", VALUE | TASKCOLOR,
251  0);
252 }
253 
254 // wrapper around QUARK_Insert_Task .. makes code more readable at task insertion time
255 void QUARK_Insert_Task_ztrmm(Quark *quark, Quark_Task_Flags *task_flags,
256  magma_int_t m,
257  magma_int_t n,
258  cuDoubleComplex alpha,
259  cuDoubleComplex *a,
260  magma_int_t lda,
261  cuDoubleComplex **b,
262  magma_int_t ldb,
263  cuDoubleComplex beta,
264  cuDoubleComplex *c,
265  magma_int_t ldc,
266  char *dag_label,
267  magma_int_t priority)
268 {
269 
270  QUARK_Insert_Task(quark, SCHED_ztrmm, task_flags,
271  sizeof(magma_int_t), &m, VALUE,
272  sizeof(magma_int_t), &n, VALUE,
273  sizeof(cuDoubleComplex), &alpha, VALUE,
274  sizeof(cuDoubleComplex)*ldb*ldb, a, INPUT,
275  sizeof(magma_int_t), &lda, VALUE,
276  sizeof(cuDoubleComplex*), b, INPUT,
277  sizeof(magma_int_t), &ldb, VALUE,
278  sizeof(cuDoubleComplex), &beta, VALUE,
279  sizeof(cuDoubleComplex)*ldb*ldb, c, INOUT | LOCALITY,
280  sizeof(magma_int_t), &ldc, VALUE,
281  sizeof(cuDoubleComplex)*ldb*ldb, NULL, SCRATCH,
282  sizeof(magma_int_t), &priority, VALUE | TASK_PRIORITY,
283  strlen(dag_label)+1, dag_label, VALUE | TASKLABEL,
284  6, "orange", VALUE | TASKCOLOR,
285  0);
286 }
287 
288 // wrapper around QUARK_Insert_Task .. makes code more readable at task insertion time
289 void QUARK_Insert_Task_zgeqrt(Quark *quark, Quark_Task_Flags *task_flags,
290  magma_int_t m,
291  magma_int_t n,
292  cuDoubleComplex *a,
293  magma_int_t lda,
294  cuDoubleComplex *t,
295  magma_int_t ldt,
296  cuDoubleComplex *tau,
297  char *dag_label)
298 {
299 
300  magma_int_t priority = 1000;
301 
302  QUARK_Insert_Task(quark, SCHED_zgeqrt, task_flags,
303  sizeof(magma_int_t), &m, VALUE,
304  sizeof(magma_int_t), &n, VALUE,
305  sizeof(magma_int_t), &ldt, VALUE,
306  sizeof(cuDoubleComplex)*m*n, a, INOUT | LOCALITY,
307  sizeof(magma_int_t), &lda, VALUE,
308  sizeof(cuDoubleComplex)*ldt*ldt, t, OUTPUT,
309  sizeof(magma_int_t), &ldt, VALUE,
310  sizeof(cuDoubleComplex)*ldt, tau, OUTPUT,
311  sizeof(cuDoubleComplex)*ldt*ldt, NULL, SCRATCH,
312  sizeof(magma_int_t), &priority, VALUE | TASK_PRIORITY,
313  strlen(dag_label)+1, dag_label, VALUE | TASKLABEL,
314  6, "green", VALUE | TASKCOLOR,
315  0);
316 
317 }
318 
319 // wrapper around QUARK_Insert_Task .. makes code more readable at task insertion time
320 void QUARK_Insert_Task_zlarfb(Quark *quark, Quark_Task_Flags *task_flags,
321  magma_int_t m,
322  magma_int_t n,
323  magma_int_t mm,
324  magma_int_t nn,
325  magma_int_t ib,
326  cuDoubleComplex *v,
327  magma_int_t ldv,
328  cuDoubleComplex *c,
329  magma_int_t ldc,
330  cuDoubleComplex *t,
331  magma_int_t ldt,
332  cuDoubleComplex **w,
333  magma_int_t ldw,
334  char *dag_label,
335  magma_int_t priority)
336 
337 {
338 
339  QUARK_Insert_Task(quark, SCHED_zlarfb, task_flags,
340  sizeof(magma_int_t), &m, VALUE,
341  sizeof(magma_int_t), &n, VALUE,
342  sizeof(magma_int_t), &mm, VALUE,
343  sizeof(magma_int_t), &nn, VALUE,
344  sizeof(magma_int_t), &ib, VALUE,
345  sizeof(cuDoubleComplex)*m*n, v, INPUT,
346  sizeof(magma_int_t), &ldv, VALUE,
347  sizeof(cuDoubleComplex)*m*n, c, INPUT,
348  sizeof(magma_int_t), &ldc, VALUE,
349  sizeof(cuDoubleComplex)*ib*ib, t, INPUT,
350  sizeof(magma_int_t), &ldt, VALUE,
351  sizeof(cuDoubleComplex*), w, OUTPUT | LOCALITY,
352  sizeof(magma_int_t), &ldw, VALUE,
353  sizeof(magma_int_t), &priority, VALUE | TASK_PRIORITY,
354  strlen(dag_label)+1, dag_label, VALUE | TASKLABEL,
355  6, "cyan", VALUE | TASKCOLOR,
356  0);
357 
358 }
359 
360 extern "C" magma_int_t
361 magma_zgeqrf_mc( magma_context *cntxt, magma_int_t *m, magma_int_t *n,
362  cuDoubleComplex *a, magma_int_t *lda, cuDoubleComplex *tau,
363  cuDoubleComplex *work, magma_int_t *lwork,
364  magma_int_t *info)
365 {
366 /* -- MAGMA (version 1.2.0) --
367  Univ. of Tennessee, Knoxville
368  Univ. of California, Berkeley
369  Univ. of Colorado, Denver
370  May 2012
371 
372  Purpose
373  =======
374 
375  ZGEQRF computes a QR factorization of a complex M-by-N matrix A:
376  A = Q * R.
377 
378  Arguments
379  =========
380  CNTXT (input) MAGMA_CONTEXT
381  CNTXT specifies the MAGMA hardware context for this routine.
382 
383  M (input) magma_int_tEGER
384  The number of rows of the matrix A. M >= 0.
385 
386  N (input) magma_int_tEGER
387  The number of columns of the matrix A. N >= 0.
388 
389  A (input/output) COMPLEX_16 array, dimension (LDA,N)
390  On entry, the M-by-N matrix A.
391  On exit, the elements on and above the diagonal of the array
392  contain the min(M,N)-by-N upper trapezoidal matrix R (R is
393  upper triangular if m >= n); the elements below the diagonal,
394  with the array TAU, represent the orthogonal matrix Q as a
395  product of min(m,n) elementary reflectors (see Further
396  Details).
397 
398  LDA (input) magma_int_tEGER
399  The leading dimension of the array A. LDA >= max(1,M).
400 
401  TAU (output) COMPLEX_16 array, dimension (min(M,N))
402  The scalar factors of the elementary reflectors (see Further
403  Details).
404 
405  WORK (workspace/output) COMPLEX_16 array, dimension (MAX(1,LWORK))
406  On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
407 
408  LWORK (input) magma_int_tEGER
409  The dimension of the array WORK. LWORK >= N*NB.
410 
411  If LWORK = -1, then a workspace query is assumed; the routine
412  only calculates the optimal size of the WORK array, returns
413  this value as the first entry of the WORK array, and no error
414  message related to LWORK is issued.
415 
416  INFO (output) magma_int_tEGER
417  = 0: successful exit
418  < 0: if INFO = -i, the i-th argument had an illegal value
419 
420  Further Details
421  ===============
422  The matrix Q is represented as a product of elementary reflectors
423 
424  Q = H(1) H(2) . . . H(k), where k = min(m,n).
425 
426  Each H(i) has the form
427 
428  H(i) = I - tau * v * v'
429 
430  where tau is a complex scalar, and v is a complex vector with
431  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
432  and tau in TAU(i).
433  ==================================================================== */
434 
435  if (cntxt->num_cores == 1 && cntxt->num_gpus == 1)
436  {
437  //magma_int_t result = magma_zgeqrf(*m, *n, a, *lda, tau, work, *lwork, info);
438  //return result;
439  }
440 
441  magma_int_t i,j,l;
442 
443  magma_int_t ii=-1,jj=-1,ll=-1;
444 
445  Quark* quark = cntxt->quark;
446 
447  // DAG labels
448  char sgeqrt_dag_label[1000];
449  char slarfb_dag_label[1000];
450  char strmm_dag_label[1000];
451  char sgemm_dag_label[1000];
452 
453  *info = 0;
454 
455  cuDoubleComplex c_one = MAGMA_Z_ONE;
456  cuDoubleComplex c_neg_one = MAGMA_Z_NEG_ONE;
457 
458  magma_int_t nb = (cntxt->nb ==-1)? magma_get_zpotrf_nb(*n): cntxt->nb;
459 
460  magma_int_t lwkopt = *n * nb;
461  work[0] = MAGMA_Z_MAKE( (double)lwkopt, 0 );
462 
463  long int lquery = *lwork == -1;
464 
465  // check input arguments
466  if (*m < 0) {
467  *info = -1;
468  } else if (*n < 0) {
469  *info = -2;
470  } else if (*lda < max(1,*m)) {
471  *info = -4;
472  } else if (*lwork < max(1,*n) && ! lquery) {
473  *info = -7;
474  }
475  if (*info != 0) {
476  magma_xerbla( __func__, -(*info) );
477  return MAGMA_ERR_ILLEGAL_VALUE;
478  }
479  else if (lquery)
480  return 0;
481 
482  magma_int_t k = min(*m,*n);
483  if (k == 0) {
484  work[0] = c_one;
485  return 0;
486  }
487 
488  magma_int_t nt = (((*n)%nb) == 0) ? (*n)/nb : (*n)/nb + 1;
489  magma_int_t mt = (((*m)%nb) == 0) ? (*m)/nb : (*m)/nb + 1;
490 
491  cuDoubleComplex **local_work = (cuDoubleComplex**) malloc(sizeof(cuDoubleComplex*)*(nt-1)*mt);
492  memset(local_work, 0, sizeof(cuDoubleComplex*)*(nt-1)*mt);
493 
494  magma_int_t priority;
495 
496  // traverse diagonal blocks
497  for (i = 0; i < k; i += nb) {
498 
499  ii++;
500 
501  jj = ii;
502 
503  sprintf(sgeqrt_dag_label, "GEQRT %d",ii);
504 
505  // factor diagonal block, also compute T matrix
506  QUARK_Insert_Task_zgeqrt(quark,
507  0, (*m)-i, min(nb,(*n)-i), A(i,i), *lda, T(i), nb, &tau[i], sgeqrt_dag_label);
508 
509  if (i > 0) {
510 
511  priority = 100;
512 
513  // update panels in a left looking fashion
514  for (j = (i-nb) + (2*nb); j < *n; j += nb) {
515 
516  jj++;
517 
518  ll = ii-1;
519 
520  sprintf(slarfb_dag_label, "LARFB %d %d",ii-1, jj);
521 
522  // perform part of update
523  QUARK_Insert_Task_zlarfb(quark, 0,
524  (*m)-(i-nb), min(nb,(*n)-(i-nb)), min(nb,(*m)-(i-nb)), min(nb,(*n)-j), nb,
525  A(i-nb,i-nb), *lda, A(i-nb,j), *lda, T(i-nb), nb, W(ii-1,jj), nb, slarfb_dag_label, priority);
526 
527  sprintf(strmm_dag_label, "TRMM %d %d",ii-1, jj);
528 
529  // perform more of update
530  QUARK_Insert_Task_ztrmm(quark, 0, min(nb,(*m)-(i-nb)), min(nb,(*n)-j), c_neg_one,
531  A(i-nb,i-nb), *lda, W(ii-1,jj), nb, c_one, A(i-nb,j), *lda, strmm_dag_label, priority);
532 
533  sprintf(sgemm_dag_label, "GEMM %d %d %d",ii-1, jj, ll);
534 
535  // finish update
536  QUARK_Insert_Task_zgemm(quark, 0, (*m)-i, min(nb,(*n)-j), min(nb,(*n)-(i-nb)), c_neg_one,
537  A(i,i-nb), *lda, W(ii-1,jj), nb, c_one, A(i,j), *lda, A(i,j), sgemm_dag_label, priority, jj);
538 
539  }
540 
541  }
542 
543  j = i + nb;
544 
545  jj = ii;
546 
547  // handle case of short wide rectangular matrix
548  if (j < (*n)) {
549 
550  priority = 0;
551 
552  jj++;
553 
554  ll = ii;
555 
556  sprintf(slarfb_dag_label, "LARFB %d %d",ii, jj);
557 
558  // perform part of update
559  QUARK_Insert_Task_zlarfb(quark, 0,
560  (*m)-i, min(nb,(*n)-i), min(nb,(*m)-i), min(nb,(*n)-j), nb,
561  A(i,i), *lda, A(i,j), *lda, T(i), nb, W(ii,jj), nb, slarfb_dag_label, priority);
562 
563  sprintf(strmm_dag_label, "TRMM %d %d",ii, jj);
564 
565  // perform more of update
566  QUARK_Insert_Task_ztrmm(quark, 0, min(nb,(*m)-i), min(nb,(*n)-j), c_neg_one,
567  A(i,i), *lda, W(ii,jj), nb, c_one, A(i,j), *lda, strmm_dag_label, priority);
568 
569  sprintf(sgemm_dag_label, "GEMM %d %d %d",ii, jj, ll);
570 
571  // finish update
572  QUARK_Insert_Task_zgemm(quark, 0, (*m)-i-nb, min(nb,(*n)-j), min(nb,(*n)-i), c_neg_one,
573  A(i+nb,i), *lda, W(ii,jj), nb, c_one, A(i+nb,j), *lda, A(i+nb,j), sgemm_dag_label, priority, jj);
574 
575  }
576 
577  }
578 
579  // wait for all tasks to finish executing
580  QUARK_Barrier(quark);
581 
582  // free memory
583  for(k = 0 ; k < (nt-1)*mt; k++) {
584  if (local_work[k] != NULL) {
585  free(local_work[k]);
586  }
587  }
588  free(local_work);
589 
590 }
591 
592 #undef A
593 #undef T
594 #undef W
595