自学教程：C++ timer_printf函数代码示例

static GENF2xq_elltrace_Harley(GEN a6, GEN T2){  pari_sp ltop = avma;  pari_timer ti;  GEN T, sqx;  GEN x, x2, t;  long n = F2x_degree(T2), N = ((n + 1)>>1) + 2;  long ispcyc;  if (n==1) return gen_m1;  if (n==2) return F2x_degree(a6) ? gen_1 : stoi(-3);  if (n==3) return F2x_degree(a6) ? (F2xq_trace(a6,T2) ?  stoi(-3): gen_1)                                  : stoi(5);  timer_start(&ti);  sqx = mkvec2(F2xq_sqrt(polx_F2x(T2[1]),T2), T2);  if (DEBUGLEVEL>1) timer_printf(&ti,"Sqrtx");  ispcyc = zx_is_pcyc(F2x_to_Flx(T2));  T = ispcyc? F2x_to_ZX(T2): F2x_canonlift(T2, N-2);  if (DEBUGLEVEL>1) timer_printf(&ti,"Teich");  T = FpX_get_red(T, int2n(N));  if (DEBUGLEVEL>1) timer_printf(&ti,"Barrett");  x = solve_AGM_eqn(F2x_to_ZX(a6), N-2, T, sqx);  if (DEBUGLEVEL>1) timer_printf(&ti,"Lift");  x2 = ZX_Z_add_shallow(ZX_shifti(x,2), gen_1);  t = (ispcyc? Z2XQ_invnorm_pcyc: Z2XQ_invnorm)(x2, T, N);  if (DEBUGLEVEL>1) timer_printf(&ti,"Norm");  if (cmpii(sqri(t), int2n(n + 2)) > 0)    t = subii(t, int2n(N));  return gerepileuptoint(ltop, t);}

void fwrite_mife_sk(const_mmap_vtable mmap, mife_sk_t sk, char *filepath) {  uint64_t t = ggh_walltime(0);  FILE *fp = fopen(filepath, "wb");  timer_printf("Starting writing Kilian matrices.../n");  fprintf(fp, "%d/n", sk->numR);  for(int i = 0; i < sk->numR; i++) {    fprintf(fp, "%ld %ld/n", sk->R[i]->r, sk->R[i]->c);    fmpz_mat_fprint_raw(fp, sk->R[i]);    fprintf(fp, "/n");    fprintf(fp, "%ld %ld/n", sk->R_inv[i]->r, sk->R_inv[i]->c);    fmpz_mat_fprint_raw(fp, sk->R_inv[i]);    fprintf(fp, "/n");  timer_printf("/r    Progress: [%lu / %lu] %8.2fs",    i, sk->numR, ggh_seconds(ggh_walltime(t)));  }  timer_printf("/n");  timer_printf("Finished writing Kilian matrices %8.2fs/n",    ggh_seconds(ggh_walltime(t)));  mmap->sk->fwrite(sk->self, fp);  fclose(fp);}

void fread_mife_sk(const_mmap_vtable mmap, mife_sk_t sk, char *filepath) {  uint64_t t = ggh_walltime(0);  FILE *fp = fopen(filepath, "rb");  timer_printf("Starting reading Kilian matrices.../n");  CHECK(fscanf(fp, "%d/n", &sk->numR), 1);  sk->R = malloc(sk->numR * sizeof(fmpz_mat_t));  sk->R_inv = malloc(sk->numR * sizeof(fmpz_mat_t));  for(int i = 0; i < sk->numR; i++) {    unsigned long r1, c1, r2, c2;    CHECK(fscanf(fp, "%lu %lu/n", &r1, &c1), 2);    fmpz_mat_init(sk->R[i], r1, c1);    fmpz_mat_fread_raw(fp, sk->R[i]);    CHECK(fscanf(fp, "/n"), 0);    CHECK(fscanf(fp, "%lu %lu/n", &r2, &c2), 2);    fmpz_mat_init(sk->R_inv[i], r2, c2);    fmpz_mat_fread_raw(fp, sk->R_inv[i]);    CHECK(fscanf(fp, "/n"), 0);  timer_printf("/r    Progress: [%lu / %lu] %8.2fs",    i, sk->numR, ggh_seconds(ggh_walltime(t)));  }  timer_printf("/n");  timer_printf("Finished reading Kilian matrices %8.2fs/n",    ggh_seconds(ggh_walltime(t)));  sk->self = malloc(mmap->sk->size);  mmap->sk->fread(sk->self, fp);  fclose(fp);}

/* Assume a = 1 [4] */static GENZ2XQ_invnorm(GEN a, GEN T, long e){  pari_timer ti;  GEN pe = int2n(e), s;  if (degpol(a)==0)    return Fp_inv(Fp_powu(gel(a,2), get_FpX_degree(T), pe), pe);  if (DEBUGLEVEL>=3) timer_start(&ti);  s = ZpXQ_log(a, T, gen_2, e);  if (DEBUGLEVEL>=3) timer_printf(&ti,"Z2XQ_log");  s = Fp_neg(FpXQ_trace(s, T, pe), pe);  if (DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_trace");  s = modii(gel(Qp_exp(cvtop(s, gen_2, e-2)),4),pe);  if (DEBUGLEVEL>=3) timer_printf(&ti,"Qp_exp");  return s;}

//.........这里部分代码省略.........    }    /* Quick return if possible */    if (n == 0) {        return *info;    }    /* Check if matrix is very small then just call LAPACK on CPU, no need for GPU */    if (n <= 128) {        #ifdef ENABLE_DEBUG        printf("--------------------------------------------------------------/n");        printf("  warning matrix too small N=%d NB=%d, calling lapack on CPU  /n", (int) n, (int) nb);        printf("--------------------------------------------------------------/n");        #endif        lapackf77_dsygvd(&itype, jobz_, uplo_,                         &n, A, &lda, B, &ldb,                         w, work, &lwork,                         iwork, &liwork, info);        *m = n;        return *info;    }    /* Form A Cholesky factorization of B. */    magma_timer_t time=0;    timer_start( time );    magma_dpotrf_m(nrgpu, uplo, n, B, ldb, info);    if (*info != 0) {        *info = n + *info;        return *info;    }    timer_stop( time );    timer_printf( "time dpotrf_m = %6.2f/n", time );    timer_start( time );    /* Transform problem to standard eigenvalue problem and solve. */    magma_dsygst_m(nrgpu, itype, uplo, n, A, lda, B, ldb, info);    timer_stop( time );    timer_printf( "time dsygst_m = %6.2f/n", time );    timer_start( time );    magma_dsyevdx_2stage_m(nrgpu, jobz, range, uplo, n, A, lda, vl, vu, il, iu, m, w, work, lwork, iwork, liwork, info);    timer_stop( time );    timer_printf( "time dsyevdx_2stage_m = %6.2f/n", time );    if (wantz && *info == 0) {        timer_start( time );        /* Backtransform eigenvectors to the original problem. */        if (itype == 1 || itype == 2) {            /* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;               backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */            if (lower) {                trans = MagmaTrans;            } else {                trans = MagmaNoTrans;            }            magma_dtrsm_m(nrgpu, MagmaLeft, uplo, trans, MagmaNonUnit, n, *m, d_one, B, ldb, A, lda);        }        else if (itype == 3) {            /* For B*A*x=(lambda)*x;               backtransform eigenvectors: x = L*y or U'*y */

//.........这里部分代码省略.........    if (n <= 128) {        lapackf77_ssygvd( &itype, jobz_, uplo_,                          &n, A, &lda, B, &ldb,                          w, work, &lwork,                          iwork, &liwork, info );        *mout = n;        return *info;    }    if (MAGMA_SUCCESS != magma_smalloc( &dA, n*ldda ) ||        MAGMA_SUCCESS != magma_smalloc( &dB, n*lddb )) {        magma_free( dA );        magma_free( dB );        *info = MAGMA_ERR_DEVICE_ALLOC;        return *info;    }    /* Form a Cholesky factorization of B. */    magma_ssetmatrix( n, n, B, ldb, dB, lddb );    magma_ssetmatrix_async( n, n,                            A,  lda,                            dA, ldda, stream );    magma_timer_t time=0;    timer_start( time );    magma_spotrf_gpu( uplo, n, dB, lddb, info );    if (*info != 0) {        *info = n + *info;        return *info;    }    timer_stop( time );    timer_printf( "time spotrf_gpu = %6.2f/n", time );    magma_queue_sync( stream );    magma_sgetmatrix_async( n, n,                            dB, lddb,                            B,  ldb, stream );    timer_start( time );    /* Transform problem to standard eigenvalue problem and solve. */    magma_ssygst_gpu( itype, uplo, n, dA, ldda, dB, lddb, info );    timer_stop( time );    timer_printf( "time ssygst_gpu = %6.2f/n", time );    /* simple fix to be able to run bigger size.     * set dB=NULL so we know to re-allocate below     * TODO: have dwork here that will be used as dB and then passed to  ssyevd.     */    if (n > 5000) {        magma_queue_sync( stream );        magma_free( dB );  dB=NULL;    }    timer_start( time );    magma_ssyevdx_gpu( jobz, range, uplo, n, dA, ldda, vl, vu, il, iu, mout, w, A, lda,                       work, lwork, iwork, liwork, info );    timer_stop( time );    timer_printf( "time ssyevdx_gpu = %6.2f/n", time );    if (wantz && *info == 0) {        timer_start( time );        

//.........这里部分代码省略.........    /* Scale matrix to allowable range, if necessary. */    anrm = magmablas_zlanhe( MagmaMaxNorm, uplo, n, dA, ldda, dwork );    iscale = 0;    sigma  = 1;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        magmablas_zlascl( uplo, 0, 0, 1., sigma, n, n, dA, ldda, info );    }    /* Call ZHETRD to reduce Hermitian matrix to tridiagonal form. */    // zhetrd rwork: e (n)    // zstedx rwork: e (n) + llrwk (1 + 4*N + 2*N**2)  ==>  1 + 5n + 2n^2    inde   = 0;    indrwk = inde + n;    llrwk  = lrwork - indrwk;    // zhetrd work: tau (n) + llwork (n*nb)  ==>  n + n*nb    // zstedx work: tau (n) + z (n^2)    // zunmtr work: tau (n) + z (n^2) + llwrk2 (n or n*nb)  ==>  2n + n^2, or n + n*nb + n^2    indtau = 0;    indwrk = indtau + n;    //indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    //llwrk2 = lwork - indwk2;    magma_timer_t time=0;    timer_start( time );#ifdef FAST_HEMV    magma_zhetrd2_gpu( uplo, n, dA, ldda, w, &rwork[inde],                       &work[indtau], wA, ldwa, &work[indwrk], llwork,                       dC, ldwork, &iinfo );#else    magma_zhetrd_gpu ( uplo, n, dA, ldda, w, &rwork[inde],                       &work[indtau], wA, ldwa, &work[indwrk], llwork,                       &iinfo );#endif    timer_stop( time );    timer_printf( "time zhetrd_gpu = %6.2f/n", time );    /* For eigenvalues only, call DSTERF.  For eigenvectors, first call       ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call ZUNMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_dsterf( &n, w, &rwork[inde], info );        magma_dmove_eig( range, n, w, &il, &iu, vl, vu, mout );    }    else {        timer_start( time );        magma_zstedx( range, n, vl, vu, il, iu, w, &rwork[inde],                      &work[indwrk], n, &rwork[indrwk],                      llrwk, iwork, liwork, dwork, info );        timer_stop( time );        timer_printf( "time zstedx = %6.2f/n", time );        timer_start( time );        magma_dmove_eig( range, n, w, &il, &iu, vl, vu, mout );        magma_zsetmatrix( n, *mout, &work[indwrk + n * (il-1) ], n, dC, lddc );        magma_zunmtr_gpu( MagmaLeft, uplo, MagmaNoTrans, n, *mout, dA, ldda, &work[indtau],                          dC, lddc, wA, ldwa, &iinfo );        magma_zcopymatrix( n, *mout, dC, lddc, dA, ldda );        timer_stop( time );        timer_printf( "time zunmtr_gpu + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        if (*info == 0) {            imax = n;        } else {            imax = *info - 1;        }        d__1 = 1. / sigma;        blasf77_dscal( &imax, &d__1, w, &ione );    }    work[0]  = MAGMA_Z_MAKE( lwmin * one_eps, 0 );  // round up    rwork[0] = lrwmin * one_eps;    iwork[0] = liwmin;    magma_queue_destroy( stream );    magma_free( dwork );    return *info;} /* magma_zheevdx_gpu */

//.........这里部分代码省略.........                    ii   = offset+jj;                    rows = maxm - ii;                    nbi  = min( nb, NBk-jj );                    for( d=0; d < ngpu; d++ ) {                        magma_setdevice(d);                                                /* wait for a block-column on GPU */                        magma_queue_sync( stream[d][0] );                                                /* start sending next column */                        if ( jj+nb < NBk ) {                            magma_zsetmatrix_async( (M-ii-nb), min(nb,NBk-jj-nb),                                                    A(ii+nb,ii+nb), lda,                                                    dA[d],          (rows-nb), stream[d][0] );                                                        /* make sure the previous update finished */                            magmablasSetKernelStream(stream[d][0]);                            //magma_queue_sync( stream[d][1] );                            magma_queue_wait_event( stream[d][0], event[d][(1+jj/nb)%2] );                                                        /* transpose next column */                            magmablas_ztranspose( M-ii-nb, nb, dA[d], rows-nb, dPT(d,0,(1+jj/nb)%2), nb );                        }                                                /* update with the block column */                        magmablasSetKernelStream(stream[d][1]);                        magma_ztrsm( MagmaRight, MagmaUpper, MagmaNoTrans, MagmaUnit,                                     n_local[d], nbi, c_one, dPT(d,0,(jj/nb)%2), nb, dAT(d,ib,0), ldn_local );                        if ( M > ii+nb ) {                            magma_zgemm( MagmaNoTrans, MagmaNoTrans,                                n_local[d], M-(ii+nb), nbi, c_neg_one, dAT(d,ib,0), ldn_local,                                dPT(d,1,(jj/nb)%2), nb, c_one, dAT(d,ib+1,0), ldn_local );                        }                        magma_event_record( event[d][(jj/nb)%2], stream[d][1] );                                        } /* end of for each block-columns in a big-panel */                }            } /* end of for each previous big-panels */            for( d=0; d < ngpu; d++ ) {                magma_setdevice(d);                magma_queue_sync( stream[d][0] );                magma_queue_sync( stream[d][1] );            magmablasSetKernelStream(NULL);            }                /* calling magma-gpu interface to panel-factorize the big panel */            if ( M > I ) {                magma_zgetrf2_mgpu(ngpu, M-I, N, nb, I, dAT, ldn_local, ipiv+I, dA, A(0,I), lda,                                   stream, &iinfo);                if ( iinfo < 0 ) {                    *info = iinfo;                    break;                } else if ( iinfo != 0 ) {                    *info = iinfo + I * NB;                    //break;                }                /* adjust pivots */                for( ii=I; ii < min(I+N,m); ii++ )                    ipiv[ii] += I;            }            time_comp += timer_stop( time );                /* download the current big panel to CPU */            timer_start( time );            magmablas_zgetmatrix_transpose_mgpu(ngpu, stream, dAT, ldn_local, A(0,I), lda, dA, maxm, M, N, nb);            for( d=0; d < ngpu; d++ ) {                magma_setdevice(d);                magma_queue_sync( stream[d][0] );                magma_queue_sync( stream[d][1] );            magmablasSetKernelStream(NULL);            }            time_get += timer_stop( time );        } /* end of for */            timer_stop( time_total );        flops = FLOPS_ZGETRF( m, n ) / 1e9;        timer_printf(" memory-allocation time: %e/n", time_alloc );        timer_printf(" NB=%d nb=%d/n", (int) NB, (int) nb );        timer_printf(" memcopy and transpose %e seconds/n", time_set );        timer_printf(" total time %e seconds/n", time_total );        timer_printf(" Performance %f GFlop/s, %f seconds without htod and dtoh/n",     flops / (time_comp),               time_comp               );        timer_printf(" Performance %f GFlop/s, %f seconds with    htod/n",              flops / (time_comp + time_set),    time_comp + time_set    );        timer_printf(" Performance %f GFlop/s, %f seconds with    dtoh/n",              flops / (time_comp + time_get),    time_comp + time_get    );        timer_printf(" Performance %f GFlop/s, %f seconds without memory-allocation/n", flops / (time_total - time_alloc), time_total - time_alloc );            for( d=0; d < ngpu0; d++ ) {            magma_setdevice(d);            magma_free( dA[d] );            magma_event_destroy( event[d][0] );            magma_event_destroy( event[d][1] );            magma_queue_destroy( stream[d][0] );            magma_queue_destroy( stream[d][1] );        }        magma_setdevice( orig_dev );        magmablasSetKernelStream( orig_stream );    }    if ( *info >= 0 )        magma_zgetrf_piv(m, n, NB, A, lda, ipiv, info);    return *info;} /* magma_zgetrf_m */

//.........这里部分代码省略.........    if (*info != 0) {        magma_xerbla( __func__, -(*info) );        return *info;    }    else if (lquery) {        return *info;    }    /* Quick return if possible */    if (n == 0) {        return *info;    }    /* If matrix is very small, then just call LAPACK on CPU, no need for GPU */    if (n <= 128) {        lapackf77_ssygvd( &itype, jobz_, uplo_,                          &n, A, &lda, B, &ldb,                          w, work, &lwork,                          iwork, &liwork, info );        return *info;    }    magma_timer_t time=0;    timer_start( time );    magma_spotrf_m( ngpu, uplo, n, B, ldb, info );    if (*info != 0) {        *info = n + *info;        return *info;    }    timer_stop( time );    timer_printf( "time spotrf = %6.2f/n", time );    timer_start( time );    /* Transform problem to standard eigenvalue problem and solve. */    magma_ssygst_m( ngpu, itype, uplo, n, A, lda, B, ldb, info );    timer_stop( time );    timer_printf( "time ssygst = %6.2f/n", time );    timer_start( time );    magma_ssyevd_m( ngpu, jobz, uplo, n, A, lda, w, work, lwork, iwork, liwork, info );    timer_stop( time );    timer_printf( "time ssyevd = %6.2f/n", time );    if (wantz && *info == 0) {        timer_start( time );        /* Backtransform eigenvectors to the original problem. */        if (itype == 1 || itype == 2) {            /* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;               backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */            if (lower) {                trans = MagmaTrans;            } else {                trans = MagmaNoTrans;            }            magma_strsm_m( ngpu, MagmaLeft, uplo, trans, MagmaNonUnit,                           n, n, d_one, B, ldb, A, lda );        }        else if (itype == 3) {            /* For B*A*x=(lambda)*x;

//.........这里部分代码省略.........    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_dsqrt(smlnum);    rmax = magma_dsqrt(bignum);    /* Scale matrix to allowable range, if necessary. */    anrm = lapackf77_zlanhe("M", uplo_, &n, A, &lda, rwork);    iscale = 0;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        lapackf77_zlascl(uplo_, &izero, &izero, &d_one, &sigma, &n, &n, A,                         &lda, info);    }    /* Call ZHETRD to reduce Hermitian matrix to tridiagonal form. */    // zhetrd rwork: e (n)    // zstedx rwork: e (n) + llrwk (1 + 4*N + 2*N**2)  ==>  1 + 5n + 2n^2    inde   = 0;    indrwk = inde + n;    llrwk  = lrwork - indrwk;    // zhetrd work: tau (n) + llwork (n*nb)  ==>  n + n*nb    // zstedx work: tau (n) + z (n^2)    // zunmtr work: tau (n) + z (n^2) + llwrk2 (n or n*nb)  ==>  2n + n^2, or n + n*nb + n^2    indtau = 0;    indwrk = indtau + n;    indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    magma_timer_t time=0;    timer_start( time );    magma_zhetrd(uplo, n, A, lda, w, &rwork[inde],                 &work[indtau], &work[indwrk], llwork, &iinfo);    timer_stop( time );    timer_printf( "time zhetrd = %6.2f/n", time );    /* For eigenvalues only, call DSTERF.  For eigenvectors, first call     ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the     tridiagonal matrix, then call ZUNMTR to multiply it to the Householder     transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_dsterf(&n, w, &rwork[inde], info);        magma_dmove_eig(range, n, w, &il, &iu, vl, vu, m);    }    else {        timer_start( time );        if (MAGMA_SUCCESS != magma_dmalloc( &dwork, 3*n*(n/2 + 1) )) {            *info = MAGMA_ERR_DEVICE_ALLOC;            return *info;        }        magma_zstedx(range, n, vl, vu, il, iu, w, &rwork[inde],                     &work[indwrk], n, &rwork[indrwk],                     llrwk, iwork, liwork, dwork, info);        magma_free( dwork );        timer_stop( time );        timer_printf( "time zstedx = %6.2f/n", time );        timer_start( time );        magma_dmove_eig(range, n, w, &il, &iu, vl, vu, m);        magma_zunmtr(MagmaLeft, uplo, MagmaNoTrans, n, *m, A, lda, &work[indtau],                     &work[indwrk + n * (il-1) ], n, &work[indwk2], llwrk2, &iinfo);        lapackf77_zlacpy("A", &n, m, &work[indwrk + n * (il-1)], &n, A, &lda);        timer_stop( time );        timer_printf( "time zunmtr + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        if (*info == 0) {            imax = n;        } else {            imax = *info - 1;        }        d__1 = 1. / sigma;        blasf77_dscal(&imax, &d__1, w, &ione);    }    work[0]  = MAGMA_Z_MAKE( lwmin * one_eps, 0.);  // round up    rwork[0] = lrwmin * one_eps;    iwork[0] = liwmin;    return *info;} /* magma_zheevdx */

//.........这里部分代码省略.........        *info = -19;    } else if (liwork < liwmin && ! lquery) {        *info = -21;    }        if (*info != 0) {        magma_xerbla( __func__, -(*info));        return *info;    }    else if (lquery) {        return *info;    }        /*  Quick return if possible */    if (n == 0) {        return *info;    }    /* Check if matrix is very small then just call LAPACK on CPU, no need for GPU */    if (n <= 128) {        #ifdef ENABLE_DEBUG        printf("--------------------------------------------------------------/n");        printf("  warning matrix too small N=%d NB=%d, calling lapack on CPU  /n", (int) n, (int) nb);        printf("--------------------------------------------------------------/n");        #endif        lapackf77_chegvd(&itype, jobz_, uplo_,                         &n, A, &lda, B, &ldb,                         w, work, &lwork,                         #if defined(PRECISION_z) || defined(PRECISION_c)                         rwork, &lrwork,                         #endif                         iwork, &liwork, info);        *m = n;        return *info;    }    magma_timer_t time=0;    timer_start( time );    magma_cpotrf_m(nrgpu, uplo, n, B, ldb, info);    if (*info != 0) {        *info = n + *info;        return *info;    }    timer_stop( time );    timer_printf("time cpotrf = %6.2f/n", time );    timer_start( time );    /* Transform problem to standard eigenvalue problem and solve. */    magma_chegst_m(nrgpu, itype, uplo, n, A, lda, B, ldb, info);    timer_stop( time );    timer_printf( "time chegst = %6.2f/n", time );    timer_start( time );    magma_cheevdx_m(nrgpu, jobz, range, uplo, n, A, lda, vl, vu, il, iu, m, w, work, lwork, rwork, lrwork, iwork, liwork, info);    timer_stop( time );    timer_printf( "time cheevd = %6.2f/n", time );    if (wantz && *info == 0) {        timer_start( time );        /* Backtransform eigenvectors to the original problem. */        if (itype == 1 || itype == 2) {            /* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;               backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */            if (lower) {                trans = MagmaConjTrans;            } else {                trans = MagmaNoTrans;            }            magma_ctrsm_m(nrgpu, MagmaLeft, uplo, trans, MagmaNonUnit,                          n, *m, c_one, B, ldb, A, lda);        }        else if (itype == 3) {            /* For B*A*x=(lambda)*x;               backtransform eigenvectors: x = L*y or U'*y */            if (lower) {                trans = MagmaNoTrans;            } else {                trans = MagmaConjTrans;            }            //magma_ctrmm(MagmaLeft, uplo, trans, MagmaNonUnit,            //            n, n, c_one, db, lddb, da, ldda);        }        timer_stop( time );        timer_printf( "time setmatrices trsm/mm + getmatrices = %6.2f/n", time );    }    work[0]  = MAGMA_C_MAKE( lwmin * one_eps, 0.);  // round up    rwork[0] = lrwmin * one_eps;    iwork[0] = liwmin;    return *info;} /* magma_chegvd_m */

//.........这里部分代码省略.........        *info = MAGMA_ERR_DEVICE_ALLOC;        return *info;    }    /* Get machine constants. */    safmin = lapackf77_slamch("Safe minimum");    eps    = lapackf77_slamch("Precision");    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_ssqrt( smlnum );    rmax = magma_ssqrt( bignum );    /* Scale matrix to allowable range, if necessary. */    anrm = magmablas_slansy( MagmaMaxNorm, uplo, n, dA, ldda, dwork, ldwork, queue );    iscale = 0;    sigma  = 1;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        magmablas_slascl( uplo, 0, 0, 1., sigma, n, n, dA, ldda, queue, info );    }    /* Call SSYTRD to reduce symmetric matrix to tridiagonal form. */    // ssytrd work: e (n) + tau (n) + llwork (n*nb)  ==>  2n + n*nb    // sstedx work: e (n) + tau (n) + z (n*n) + llwrk2 (1 + 4*n + n^2)  ==>  1 + 6n + 2n^2    inde   = 0;    indtau = inde   + n;    indwrk = indtau + n;    indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    magma_timer_t time=0;    timer_start( time );#ifdef FAST_SYMV    magma_ssytrd2_gpu( uplo, n, dA, ldda, w, &work[inde],                       &work[indtau], wA, ldwa, &work[indwrk], llwork,                       dwork, ldwork, &iinfo );#else    magma_ssytrd_gpu(  uplo, n, dA, ldda, w, &work[inde],                       &work[indtau], wA, ldwa, &work[indwrk], llwork,                       &iinfo );#endif    timer_stop( time );    #ifdef FAST_SYMV    timer_printf( "time ssytrd2 = %6.2f/n", time );    #else    timer_printf( "time ssytrd = %6.2f/n", time );    #endif    /* For eigenvalues only, call SSTERF.  For eigenvectors, first call       SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call SORMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_ssterf( &n, w, &work[inde], info );    }    else {        timer_start( time );        magma_sstedx( MagmaRangeAll, n, 0., 0., 0, 0, w, &work[inde],                      &work[indwrk], n, &work[indwk2],                      llwrk2, iwork, liwork, dwork, info );        timer_stop( time );        timer_printf( "time sstedx = %6.2f/n", time );        timer_start( time );        magma_ssetmatrix( n, n, &work[indwrk], n, dwork, lddc, queue );        magma_sormtr_gpu( MagmaLeft, uplo, MagmaNoTrans, n, n, dA, ldda, &work[indtau],                          dwork, lddc, wA, ldwa, &iinfo );        magma_scopymatrix( n, n, dwork, lddc, dA, ldda, queue );        timer_stop( time );        timer_printf( "time sormtr + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        d__1 = 1. / sigma;        blasf77_sscal( &n, &d__1, w, &ione );    }    work[0]  = magma_smake_lwork( lwmin );    iwork[0] = liwmin;    magma_queue_destroy( queue );    magma_free( dwork );    return *info;} /* magma_ssyevd_gpu */

//.........这里部分代码省略.........                    magma_sgemm( MagmaNoTrans, MagmaConjTrans,                                 jb, jj, nb,                                 c_neg_one, dT(d, J+jj, 0), ldda,                                            dT(d, J,    0), ldda,                                 c_one,     dA(d, J2,   0), lddla);                                        magma_ssyrk(MagmaLower, MagmaNoTrans, jb, nb,                                d_neg_one, dT(d, J+jj, 0), ldda,                                d_one,     dA(d, J2,  jj), lddla);                }                                if ( n > J+JB ) { /* off-diagonal */                    for( d=0; d < num_gpus; d++ ) {                        magma_setdevice(d);                        magmablasSetKernelStream(stream[d][2]);                                                /* local number of columns in the big panel */                        n_local[d] = (((n-J)/nb)/num_gpus)*nb;                        if (d < ((n-J)/nb)%num_gpus)                            n_local[d] += nb;                        else if (d == ((n-J)/nb)%num_gpus)                            n_local[d] += (n-J)%nb;                                                /* subtracting local number of columns in diagonal */                        J2 = nb*(JB/(nb*num_gpus));                        if ( d < (JB/nb)%num_gpus )                            J2 += nb;                        n_local[d] -= J2;                                                magma_sgemm( MagmaNoTrans, MagmaConjTrans,                                     n_local[d], JB, nb,                                     c_neg_one, dT(d, J+JB, 0), ldda,                                                dT(d, J,    0), ldda,                                     c_one,     dA(d, J2,   0), lddla);                    }                }                /* wait for the previous updates */                for( d=0; d < num_gpus; d++ ) {                    magma_setdevice(d);                    for( jj=0; jj < 3; jj++ )                        magma_queue_sync( stream[d][jj] );                    magmablasSetKernelStream(NULL);                }                magma_setdevice(0);            }                        /* factor the big panel */            h = (JB+nb-1)/nb; // big diagonal of big panel will be on CPU            // using two streams            //magma_spotrf2_mgpu(num_gpus, uplo, n-J, JB, J, J, nb,            //                   dwork, lddla, dt, ldda, A, lda, h, stream, event, &iinfo);            // using three streams            magma_spotrf3_mgpu(num_gpus, uplo, n-J, JB, J, J, nb,                               dwork, lddla, dt, ldda, A, lda, h, stream, event, &iinfo);            if ( iinfo != 0 ) {                *info = J+iinfo;                break;            }            time_sum += timer_stop( time );                        /* upload the off-diagonal big panel */            magma_sdtohpo( num_gpus, uplo, n, JB, J, J, nb, JB, A, lda, dwork, lddla, stream, &iinfo);                } /* end of for J */    } /* if upper */    } /* if nb */    timer_stop( time_total );        if ( num_gpus0 > n/nb ) {        num_gpus = n/nb;        if ( n%nb != 0 ) num_gpus ++;    } else {        num_gpus = num_gpus0;    }    for (d=0; d < num_gpus; d++ ) {        magma_setdevice(d);        for( j=0; j < 3; j++ ) {            magma_queue_destroy( stream[d][j] );        }        magma_free( dt[d] );        for( j=0; j < 5; j++ ) {            magma_event_destroy( event[d][j] );        }    }    magma_setdevice( orig_dev );    magmablasSetKernelStream( orig_stream );                     timer_printf( "/n n=%d NB=%d nb=%d/n", (int) n, (int) NB, (int) nb );    timer_printf( " Without memory allocation: %f / %f = %f GFlop/s/n",                  FLOPS_SPOTRF(n) / 1e9,  time_total,                  FLOPS_SPOTRF(n) / 1e9 / time_total );    timer_printf( " Performance %f / %f = %f GFlop/s/n",                  FLOPS_SPOTRF(n) / 1e9,  time_sum,                  FLOPS_SPOTRF(n) / 1e9 / time_sum );        return *info;} /* magma_spotrf_ooc */

//.........这里部分代码省略.........    /* Set the GPU number that holds the last panel */    id = s%num_gpus;        /* Set the local index where the last panel is */    i_local = s/num_gpus;        /* size of the last diagonal-block */    nb0 = min(m - s*nb, n - s*nb);    rows = m    - s*nb;    cols = maxm - s*nb;        if ( nb0 > 0 ) {        magma_setdevice(id);                /* wait for the last panel on cpu */        magma_queue_sync( streaml[id][1] );            /* factor on cpu */        lapackf77_sgetrf( &rows, &nb0, W(s), &ldw, ipiv+s*nb, &iinfo);        if ( (*info == 0) && (iinfo > 0) )            *info = iinfo + s*nb;                /* send the factor to gpus */        for( d=0; d < num_gpus; d++ ) {            magma_setdevice(d);            i_local2 = i_local;            if ( d < id ) i_local2 ++;                        if ( d == id || n_local[d] > i_local2*nb ) {                magma_ssetmatrix_async( rows, nb0,                                        W(s),     ldw,                                        &d_lAP[d][(s%h)*nb*maxm],                                        cols, streaml[d][1] );            }        }                for( d=0; d < num_gpus; d++ ) {            magma_setdevice(d);            magmablasSetKernelStream(streaml[d][0]);            if ( d == 0 )                magmablas_spermute_long2( lddat, dAT(d,0,0), lddat, ipiv, nb0, s*nb );            else                magmablas_spermute_long3(        dAT(d,0,0), lddat, ipiv, nb0, s*nb );        }                for( d=0; d < num_gpus; d++ ) {            magma_setdevice(d);            magmablasSetKernelStream(streaml[d][1]);                        /* wait for the pivoting to be done */            magma_queue_sync( streaml[d][0] );                        i_local2 = i_local;            if ( d < id ) i_local2++;            if ( d == id ) {                /* the panel belond to this gpu */                panel_local[d] = dAT(d,s,i_local);                                /* next column */                nb1 = n_local[d] - i_local*nb-nb0;                                magmablas_stranspose( rows, nb0, &d_lAP[d][(s%h)*nb*maxm], cols, panel_local[d], lddat );                                if ( nb1 > 0 ) {                    magma_strsm( MagmaRight, MagmaUpper, MagmaNoTrans, MagmaUnit,                                 nb1, nb0, c_one,                                 panel_local[d],        lddat,                                 dAT(d,s,i_local)+nb0, lddat);                }            } else if ( n_local[d] > i_local2*nb ) {                /* the panel belong to another gpu */                panel_local[d] = d_panel[d];                                /* next column */                nb1 = n_local[d] - i_local2*nb;                                magmablas_stranspose( rows, nb0, &d_lAP[d][(s%h)*nb*maxm], cols, panel_local[d], nb );                magma_strsm( MagmaRight, MagmaUpper, MagmaNoTrans, MagmaUnit,                             nb1, nb0, c_one,                             panel_local[d],     nb,                             dAT(d,s,i_local2), lddat);            }        }    } /* if ( nb0 > 0 ) */        /* clean up */    trace_finalize( "sgetrf_mgpu.svg","trace.css" );    for( d=0; d < num_gpus; d++ ) {        magma_setdevice(d);        magma_queue_sync( streaml[d][0] );        magma_queue_sync( streaml[d][1] );        magmablasSetKernelStream(NULL);    }    magma_setdevice(0);        timer_start( time );    timer_printf("/n Performance %f GFlop/s/n", FLOPS_SGETRF(m,n) / 1e9 / time );    return *info;} /* magma_sgetrf2_mgpu */

//.........这里部分代码省略.........                    blasf77_zgemm( "n", "n", &n, &nb2, &n2, &c_one,                                   VR, &ldvr,                                   work(0,iv   ), &n, &c_zero,                                   work(0,nb+iv), &n );                    time_gemm_sum += timer_stop( time_gemm );                                        // normalize vectors                    // TODO if somev, should copy vectors individually to correct location.                    for( k = iv; k <= nb; ++k ) {                        ii = blasf77_izamax( &n, work(0,nb+k), &ione ) - 1;                        remax = 1. / MAGMA_Z_ABS1( *work(ii,nb+k) );                        blasf77_zdscal( &n, &remax, work(0,nb+k), &ione );                    }                    lapackf77_zlacpy( "F", &n, &nb2, work(0,nb+iv), &n, VR(0,ki), &ldvr );                    iv = nb;                    timer_start( time_trsv );                }                else {                    iv -= 1;                }            } // blocked back-transform            // Restore the original diagonal elements of T.            for( k=0; k <= ki - 1; ++k ) {                *T(k,k) = *work(k,0);            }            is -= 1;        }    }    timer_stop( time_trsv );    timer_stop( time_total );    timer_printf( "trevc trsv %.4f, gemm %.4f, gemv %.4f, total %.4f/n",                  time_trsv_sum, time_gemm_sum, time_gemv_sum, time_total );    if ( leftv ) {        // ============================================================        // Compute left eigenvectors.        // iv is index of column in current block.        // Non-blocked version always uses iv=1;        // blocked     version starts with iv=1, goes up to nb.        // (Note the "0-th" column is used to store the original diagonal.)        iv = 1;        is = 0;        for( ki=0; ki < n; ++ki ) {            if ( somev ) {                if ( ! select[ki] ) {                    continue;                }            }            smin = max( ulp*MAGMA_Z_ABS1( *T(ki,ki) ), smlnum );            // --------------------------------------------------------            // Complex left eigenvector            *work(ki,iv) = c_one;            // Form right-hand side.            for( k = ki + 1; k < n; ++k ) {                *work(k,iv) = -MAGMA_Z_CNJG( *T(ki,k) );            }            // Solve conjugate-transposed triangular system:            // [ T(ki+1:n,ki+1:n) - T(ki,ki) ]**H * X = scale*work.            for( k = ki + 1; k < n; ++k ) {                *T(k,k) -= *T(ki,ki);

//.........这里部分代码省略.........        #ifdef ENABLE_DEBUG        printf("--------------------------------------------------------------/n");        printf("  warning matrix too small N=%d NB=%d, calling lapack on CPU  /n", (int) n, (int) nb);        printf("--------------------------------------------------------------/n");        #endif        lapackf77_dsyevd(jobz_, uplo_,                         &n, a, &lda,                         w, work, &lwork,                         iwork, &liwork, info);        return *info;    }    /* Get machine constants. */    safmin = lapackf77_dlamch("Safe minimum");    eps    = lapackf77_dlamch("Precision");    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_dsqrt(smlnum);    rmax = magma_dsqrt(bignum);    /* Scale matrix to allowable range, if necessary. */    anrm = lapackf77_dlansy("M", uplo_, &n, a, &lda, work );    iscale = 0;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        lapackf77_dlascl(uplo_, &izero, &izero, &d_one, &sigma, &n, &n, a,                &lda, info);    }    /* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */    // dsytrd work: e (n) + tau (n) + llwork (n*nb)  ==>  2n + n*nb    // dstedx work: e (n) + tau (n) + z (n*n) + llwrk2 (1 + 4*n + n^2)  ==>  1 + 6n + 2n^2    inde   = 0;    indtau = inde   + n;    indwrk = indtau + n;    indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    magma_timer_t time;    timer_start( time );    magma_dsytrd(uplo, n, a, lda, w, &work[inde],                 &work[indtau], &work[indwrk], llwork, queue, &iinfo);        timer_stop( time );    timer_printf( "time dsytrd = %6.2f/n", time );    /* For eigenvalues only, call DSTERF.  For eigenvectors, first call       DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call DORMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_dsterf(&n, w, &work[inde], info);    }    else {        timer_start( time );        if (MAGMA_SUCCESS != magma_dmalloc( &dwork, 3*n*(n/2 + 1) )) {            *info = MAGMA_ERR_DEVICE_ALLOC;            return *info;        }                // TTT Possible bug for n < 128        magma_dstedx(MagmaRangeAll, n, 0., 0., 0, 0, w, &work[inde],                     &work[indwrk], n, &work[indwk2],                     llwrk2, iwork, liwork, dwork, queue, info);                magma_free( dwork );                timer_stop( time );        timer_printf( "time dstedx = %6.2f/n", time );        timer_start( time );                magma_dormtr(MagmaLeft, uplo, MagmaNoTrans, n, n, a, lda, &work[indtau],                     &work[indwrk], n, &work[indwk2], llwrk2, queue, &iinfo);                lapackf77_dlacpy("A", &n, &n, &work[indwrk], &n, a, &lda);        timer_stop( time );        timer_printf( "time dormtr + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        d__1 = 1. / sigma;        blasf77_dscal(&n, &d__1, w, &ione);    }    work[0]  = lwmin * one_eps;  // round up    iwork[0] = liwmin;    return *info;} /* magma_dsyevd */

//.........这里部分代码省略.........#pragma omp barrier                //reduce the number of used threads to have enough S workspace                tot = min(n1, omp_get_num_threads());                if (id < tot) {                    ib = (  id   * rk) / tot + iil - 1;                    ie = ((id+1) * rk) / tot + iil - 1;                    ik = ie - ib;                }                else {                    ib = -1;                    ie = -1;                    ik = -1;                }                // Compute eigenvectors of the modified rank-1 modification.                for (j = ib; j < ie; ++j) {                    for (i = 0; i < k; ++i)                        s[id*k + i] = w[i] / *Q(i,j);                    temp = magma_cblas_dnrm2( k, s+id*k, 1 );                    for (i = 0; i < k; ++i) {                        magma_int_t iii = indx[i] - 1;                        *Q(i,j) = s[id*k + iii] / temp;                    }                }            }        }    }    if (*info != 0)        return *info;    timer_stop( time );    timer_printf( "eigenvalues/vector D+zzT = %6.2f/n", time );#else    /////////////////////////////////////////////////////////////////////////////////    // Non openmp implementation    /////////////////////////////////////////////////////////////////////////////////    magma_timer_t time=0;    timer_start( time );    for (i = 0; i < k; ++i)        dlamda[i]=lapackf77_dlamc3(&dlamda[i], &dlamda[i]) - dlamda[i];    for (j = 0; j < k; ++j) {        magma_int_t tmpp=j+1;        magma_int_t iinfo = 0;        lapackf77_dlaed4(&k, &tmpp, dlamda, w, Q(0,j), &rho, &d[j], &iinfo);        // If the zero finder fails, the computation is terminated.        if (iinfo != 0)            *info=iinfo;    }    if (*info != 0)        return *info;    //Prepare the INDXQ sorting permutation.    magma_int_t nk = n - k;    lapackf77_dlamrg( &k, &nk, d, &ione, &ineg_one, indxq);    //compute the lower and upper bound of the non-deflated eigenvectors    if (valeig)        magma_dvrange(k, d, &iil, &iiu, vl, vu);    else if (indeig)        magma_dirange(k, indxq, &iil, &iiu, il, iu);    else {

//.........这里部分代码省略.........        lapackf77_dsyevd(jobz_, uplo_,                         &n, A, &lda,                         w, work, &lwork,                         iwork, &liwork, info);        return *info;    }    /* Get machine constants. */    safmin = lapackf77_dlamch("Safe minimum");    eps = lapackf77_dlamch("Precision");    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_dsqrt(smlnum);    rmax = magma_dsqrt(bignum);    /* Scale matrix to allowable range, if necessary. */    anrm = lapackf77_dlansy("M", uplo_, &n, A, &lda, work);    iscale = 0;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        lapackf77_dlascl(uplo_, &izero, &izero, &d_one, &sigma, &n, &n, A,                &lda, info);    }    /* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */    // dsytrd work: e (n) + tau (n) + llwork (n*nb)  ==>  2n + n*nb    // dstedx work: e (n) + tau (n) + z (n*n) + llwrk2 (1 + 4*n + n^2)  ==>  1 + 6n + 2n^2    inde   = 0;    indtau = inde   + n;    indwrk = indtau + n;    indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    magma_timer_t time=0;    timer_start( time );    magma_dsytrd_mgpu(nrgpu, 1, uplo, n, A, lda, w, &work[inde],                      &work[indtau], &work[indwrk], llwork, &iinfo);    timer_stop( time );    timer_printf( "time dsytrd = %6.2f/n", time );    /* For eigenvalues only, call DSTERF.  For eigenvectors, first call       DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call DORMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_dsterf(&n, w, &work[inde], info);    }    else {        timer_start( time );#ifdef USE_SINGLE_GPU        if (MAGMA_SUCCESS != magma_dmalloc( &dwork, 3*n*(n/2 + 1) )) {            *info = MAGMA_ERR_DEVICE_ALLOC;            return *info;        }        magma_dstedx(MagmaRangeAll, n, 0., 0., 0, 0, w, &work[inde],                     &work[indwrk], n, &work[indwk2],                     llwrk2, iwork, liwork, dwork, info);        magma_free( dwork );#else        magma_dstedx_m(nrgpu, MagmaRangeAll, n, 0., 0., 0, 0, w, &work[inde],                       &work[indwrk], n, &work[indwk2],                       llwrk2, iwork, liwork, info);#endif        timer_stop( time );        timer_printf( "time dstedc = %6.2f/n", time );        timer_start( time );        magma_dormtr_m(nrgpu, MagmaLeft, uplo, MagmaNoTrans, n, n, A, lda, &work[indtau],                       &work[indwrk], n, &work[indwk2], llwrk2, &iinfo);        lapackf77_dlacpy("A", &n, &n, &work[indwrk], &n, A, &lda);        timer_stop( time );        timer_printf( "time dormtr + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        d__1 = 1. / sigma;        blasf77_dscal(&n, &d__1, w, &ione);    }    work[0]  = lwmin * one_eps;  // round up    iwork[0] = liwmin;    return *info;} /* magma_dsyevd_m */

//.........这里部分代码省略.........    if (wantvl) {        /* Undo balancing of left eigenvectors         * (Workspace: need N) */        lapackf77_sgebak( "B", "L", &n, &ilo, &ihi, &work[ibal], &n,                          VL, &ldvl, &ierr );        /* Normalize left eigenvectors and make largest component real */        for (i = 0; i < n; ++i) {            if ( wi[i] == 0. ) {                scl = 1. / magma_cblas_snrm2( n, VL(0,i), 1 );                blasf77_sscal( &n, &scl, VL(0,i), &ione );            }            else if ( wi[i] > 0. ) {                d__1 = magma_cblas_snrm2( n, VL(0,i), 1 );                d__2 = magma_cblas_snrm2( n, VL(0,i+1), 1 );                scl = 1. / lapackf77_slapy2( &d__1, &d__2 );                blasf77_sscal( &n, &scl, VL(0,i), &ione );                blasf77_sscal( &n, &scl, VL(0,i+1), &ione );                for (k = 0; k < n; ++k) {                    /* Computing 2nd power */                    d__1 = *VL(k,i);                    d__2 = *VL(k,i+1);                    work[iwrk + k] = d__1*d__1 + d__2*d__2;                }                k = blasf77_isamax( &n, &work[iwrk], &ione ) - 1;  // subtract 1; k is 0-based                lapackf77_slartg( VL(k,i), VL(k,i+1), &cs, &sn, &r );                blasf77_srot( &n, VL(0,i), &ione, VL(0,i+1), &ione, &cs, &sn );                *VL(k,i+1) = 0.;            }        }    }    if (wantvr) {        /* Undo balancing of right eigenvectors         * (Workspace: need N) */        lapackf77_sgebak( "B", "R", &n, &ilo, &ihi, &work[ibal], &n,                          VR, &ldvr, &ierr );        /* Normalize right eigenvectors and make largest component real */        for (i = 0; i < n; ++i) {            if ( wi[i] == 0. ) {                scl = 1. / magma_cblas_snrm2( n, VR(0,i), 1 );                blasf77_sscal( &n, &scl, VR(0,i), &ione );            }            else if ( wi[i] > 0. ) {                d__1 = magma_cblas_snrm2( n, VR(0,i), 1 );                d__2 = magma_cblas_snrm2( n, VR(0,i+1), 1 );                scl = 1. / lapackf77_slapy2( &d__1, &d__2 );                blasf77_sscal( &n, &scl, VR(0,i), &ione );                blasf77_sscal( &n, &scl, VR(0,i+1), &ione );                for (k = 0; k < n; ++k) {                    /* Computing 2nd power */                    d__1 = *VR(k,i);                    d__2 = *VR(k,i+1);                    work[iwrk + k] = d__1*d__1 + d__2*d__2;                }                k = blasf77_isamax( &n, &work[iwrk], &ione ) - 1;  // subtract 1; k is 0-based                lapackf77_slartg( VR(k,i), VR(k,i+1), &cs, &sn, &r );                blasf77_srot( &n, VR(0,i), &ione, VR(0,i+1), &ione, &cs, &sn );                *VR(k,i+1) = 0.;            }        }    }CLEANUP:    /* Undo scaling if necessary */    if (scalea) {        // converged eigenvalues, stored in wr[i+1:n] and wi[i+1:n] for i = INFO        magma_int_t nval = n - (*info);        magma_int_t ld = max( nval, 1 );        lapackf77_slascl( "G", &izero, &izero, &cscale, &anrm, &nval, &ione, wr + (*info), &ld, &ierr );        lapackf77_slascl( "G", &izero, &izero, &cscale, &anrm, &nval, &ione, wi + (*info), &ld, &ierr );        if (*info > 0) {            // first ilo columns were already upper triangular,            // so the corresponding eigenvalues are also valid.            nval = ilo - 1;            lapackf77_slascl( "G", &izero, &izero, &cscale, &anrm, &nval, &ione, wr, &n, &ierr );            lapackf77_slascl( "G", &izero, &izero, &cscale, &anrm, &nval, &ione, wi, &n, &ierr );        }    }    #if defined(Version3)    magma_free( dT );    #endif    #if defined(Version5)    magma_free_cpu( T );    #endif        timer_stop( time_total );    flops_stop( flop_total );    timer_printf( "sgeev times n %5d, gehrd %7.3f, unghr %7.3f, hseqr %7.3f, trevc %7.3f, total %7.3f, sum %7.3f/n",                  (int) n, time_gehrd, time_unghr, time_hseqr, time_trevc, time_total, time_sum );    timer_printf( "sgeev flops n %5d, gehrd %7lld, unghr %7lld, hseqr %7lld, trevc %7lld, total %7lld, sum %7lld/n",                  (int) n, flop_gehrd, flop_unghr, flop_hseqr, flop_trevc, flop_total, flop_sum );        work[0] = magma_smake_lwork( optwrk );        return *info;} /* magma_sgeev */

//.........这里部分代码省略.........        return *info;    }    /* Get machine constants. */    safmin = lapackf77_slamch("Safe minimum");    eps    = lapackf77_slamch("Precision");    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_ssqrt(smlnum);    rmax = magma_ssqrt(bignum);    /* Scale matrix to allowable range, if necessary. */    anrm = magmablas_slansy(MagmaMaxNorm, uplo, n, dA, ldda, dwork);    iscale = 0;    sigma  = 1;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        magmablas_slascl(uplo, 0, 0, 1., sigma, n, n, dA, ldda, info);    }    /* Call SSYTRD to reduce symmetric matrix to tridiagonal form. */    // ssytrd work: e (n) + tau (n) + llwork (n*nb)  ==>  2n + n*nb    // sstedx work: e (n) + tau (n) + z (n*n) + llwrk2 (1 + 4*n + n^2)  ==>  1 + 6n + 2n^2    inde   = 0;    indtau = inde   + n;    indwrk = indtau + n;    indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    magma_timer_t time=0;    timer_start( time );#ifdef FAST_SYMV    magma_ssytrd2_gpu(uplo, n, dA, ldda, w, &work[inde],                      &work[indtau], wA, ldwa, &work[indwrk], llwork,                      dwork, n*lddc, &iinfo);#else    magma_ssytrd_gpu(uplo, n, dA, ldda, w, &work[inde],                     &work[indtau], wA, ldwa, &work[indwrk], llwork,                     &iinfo);#endif    timer_stop( time );    timer_printf( "time ssytrd = %6.2f/n", time );    /* For eigenvalues only, call SSTERF.  For eigenvectors, first call       SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call SORMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_ssterf(&n, w, &work[inde], info);        magma_smove_eig(range, n, w, &il, &iu, vl, vu, m);    }    else {        timer_start( time );        magma_sstedx(range, n, vl, vu, il, iu, w, &work[inde],                     &work[indwrk], n, &work[indwk2],                     llwrk2, iwork, liwork, dwork, info);        timer_stop( time );        timer_printf( "time sstedx = %6.2f/n", time );        timer_start( time );        magma_smove_eig(range, n, w, &il, &iu, vl, vu, m);        magma_ssetmatrix( n, *m, &work[indwrk + n* (il-1) ], n, dwork, lddc );        magma_sormtr_gpu(MagmaLeft, uplo, MagmaNoTrans, n, *m, dA, ldda, &work[indtau],                         dwork, lddc, wA, ldwa, &iinfo);        magma_scopymatrix( n, *m, dwork, lddc, dA, ldda );        timer_stop( time );        timer_printf( "time sormtr + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        d__1 = 1. / sigma;        blasf77_sscal(&n, &d__1, w, &ione);    }    work[0]  = lwmin * one_eps;  // round up    iwork[0] = liwmin;    magma_queue_destroy( stream );    magma_free( dwork );    return *info;} /* magma_ssyevd_gpu */

//.........这里部分代码省略.........    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_ssqrt(smlnum);    rmax = magma_ssqrt(bignum);    /* Scale matrix to allowable range, if necessary. */    anrm = magmablas_clanhe(MagmaMaxNorm, uplo, n, dA, ldda, dwork);    iscale = 0;    sigma  = 1;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        magmablas_clascl(uplo, 0, 0, 1., sigma, n, n, dA, ldda, info);    }    /* Call CHETRD to reduce Hermitian matrix to tridiagonal form. */    // chetrd rwork: e (n)    // cstedx rwork: e (n) + llrwk (1 + 4*N + 2*N**2)  ==>  1 + 5n + 2n^2    inde   = 0;    indrwk = inde + n;    llrwk  = lrwork - indrwk;    // chetrd work: tau (n) + llwork (n*nb)  ==>  n + n*nb    // cstedx work: tau (n) + z (n^2)    // cunmtr work: tau (n) + z (n^2) + llwrk2 (n or n*nb)  ==>  2n + n^2, or n + n*nb + n^2    indtau = 0;    indwrk = indtau + n;    indwk2 = indwrk + n*n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    magma_timer_t time=0;    timer_start( time );#ifdef FAST_HEMV    magma_chetrd2_gpu(uplo, n, dA, ldda, w, &rwork[inde],                      &work[indtau], wA, ldwa, &work[indwrk], llwork,                      dC, n*lddc, &iinfo);#else    magma_chetrd_gpu (uplo, n, dA, ldda, w, &rwork[inde],                      &work[indtau], wA, ldwa, &work[indwrk], llwork, &iinfo);#endif    timer_stop( time );    timer_printf( "time chetrd_gpu = %6.2f/n", time );    /* For eigenvalues only, call SSTERF.  For eigenvectors, first call       CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call CUNMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_ssterf(&n, w, &rwork[inde], info);    }    else {        timer_start( time );        magma_cstedx( MagmaRangeAll, n, 0., 0., 0, 0, w, &rwork[inde],                      &work[indwrk], n, &rwork[indrwk],                      llrwk, iwork, liwork, dwork, info);        timer_stop( time );        timer_printf( "time cstedx = %6.2f/n", time );        timer_start( time );        magma_csetmatrix( n, n, &work[indwrk], n, dC, lddc );        magma_cunmtr_gpu(MagmaLeft, uplo, MagmaNoTrans, n, n, dA, ldda, &work[indtau],                         dC, lddc, wA, ldwa, &iinfo);        magma_ccopymatrix( n, n, dC, lddc, dA, ldda );        timer_stop( time );        timer_printf( "time cunmtr_gpu + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        if (*info == 0) {            imax = n;        } else {            imax = *info - 1;        }        d__1 = 1. / sigma;        blasf77_sscal(&imax, &d__1, w, &ione);    }    work[0]  = MAGMA_C_MAKE( lwmin * one_eps, 0.);  // round up    rwork[0] = lrwmin * one_eps;    iwork[0] = liwmin;    magma_queue_destroy( stream );    magma_free( dwork );    return *info;} /* magma_cheevd_gpu */

//.........这里部分代码省略.........#pragma omp barrier                //reduce the number of used threads to have enough S workspace                tot = min(n1, omp_get_num_threads());                if (id < tot) {                    ib = (  id   * rk) / tot + iil - 1;                    ie = ((id+1) * rk) / tot + iil - 1;                    ik = ie - ib;                }                else {                    ib = -1;                    ie = -1;                    ik = -1;                }                // Compute eigenvectors of the modified rank-1 modification.                for (j = ib; j < ie; ++j) {                    for (i = 0; i < k; ++i)                        s[id*k + i] = w[i] / *Q(i,j);                    temp = magma_cblas_snrm2( k, s+id*k, 1 );                    for (i = 0; i < k; ++i) {                        magma_int_t iii = indx[i] - 1;                        *Q(i,j) = s[id*k + iii] / temp;                    }                }            }        }    }    if (*info != 0)        return *info;    timer_stop( time );    timer_printf( "eigenvalues/vector D+zzT = %6.2f/n", time );#else    /////////////////////////////////////////////////////////////////////////////////    // Non openmp implementation    /////////////////////////////////////////////////////////////////////////////////    magma_timer_t time=0;    timer_start( time );    for (i = 0; i < k; ++i)        dlamda[i]=lapackf77_slamc3(&dlamda[i], &dlamda[i]) - dlamda[i];    for (j = 0; j < k; ++j) {        magma_int_t tmpp=j+1;        magma_int_t iinfo = 0;        lapackf77_slaed4(&k, &tmpp, dlamda, w, Q(0,j), &rho, &d[j], &iinfo);        // If the zero finder fails, the computation is terminated.        if (iinfo != 0)            *info=iinfo;    }    if (*info != 0)        return *info;    //Prepare the INDXQ sorting permutation.    magma_int_t nk = n - k;    lapackf77_slamrg( &k, &nk, d, &ione, &ineg_one, indxq);    //compute the lower and upper bound of the non-deflated eigenvectors    if (valeig)        magma_svrange(k, d, &iil, &iiu, vl, vu);    else if (indeig)        magma_sirange(k, indxq, &iil, &iiu, il, iu);    else {

//.........这里部分代码省略.........                         &n, A, &lda, B, &ldb,                         w, work, &lwork,                         #if defined(PRECISION_z) || defined(PRECISION_c)                         rwork, &lrwork,                         #endif                         iwork, &liwork, info);        *m = n;        return *info;    }    // TODO: fix memory leak    if (MAGMA_SUCCESS != magma_zmalloc( &dA, n*ldda ) ||        MAGMA_SUCCESS != magma_zmalloc( &dB, n*lddb )) {        *info = MAGMA_ERR_DEVICE_ALLOC;        return *info;    }    /* Form a Cholesky factorization of B. */    magma_zsetmatrix( n, n, B, ldb, dB, lddb );    magma_zsetmatrix_async( n, n,                           A,  lda,                           dA, ldda, stream );    magma_timer_t time=0;    timer_start( time );    magma_zpotrf_gpu(uplo, n, dB, lddb, info);    if (*info != 0) {        *info = n + *info;        return *info;    }    timer_stop( time );    timer_printf( "time zpotrf_gpu = %6.2f/n", time );    magma_queue_sync( stream );    magma_zgetmatrix_async( n, n,                           dB, lddb,                           B,  ldb, stream );    /* Transform problem to standard eigenvalue problem and solve. */    timer_start( time );    magma_zhegst_gpu(itype, uplo, n, dA, ldda, dB, lddb, info);    timer_stop( time );    timer_printf( "time zhegst_gpu = %6.2f/n", time );    magma_zgetmatrix( n, n, dA, ldda, A, lda );    magma_queue_sync( stream );    magma_free( dA );    magma_free( dB );    timer_start( time );    magma_zheevdx_2stage(jobz, range, uplo, n, A, lda, vl, vu, il, iu, m, w, work, lwork, rwork, lrwork, iwork, liwork, info);    timer_stop( time );    timer_printf( "time zheevdx_2stage = %6.2f/n", time );    if (wantz && *info == 0) {        // TODO fix memory leak        if (MAGMA_SUCCESS != magma_zmalloc( &dA, n*ldda ) ||            MAGMA_SUCCESS != magma_zmalloc( &dB, n*lddb )) {            *info = MAGMA_ERR_DEVICE_ALLOC;            return *info;        }        timer_start( time );

//.........这里部分代码省略.........        lapackf77_ztrevc( lapack_side_const(side), "B", select, &n, A, &lda, VL, &ldvl,                          VR, &ldvr, &n, &nout, &work[iwrk], &rwork[irwork], &ierr );        #elif TREVC_VERSION == 2        liwrk = lwork - iwrk;        lapackf77_ztrevc3( lapack_side_const(side), "B", select, &n, A, &lda, VL, &ldvl,                           VR, &ldvr, &n, &nout, &work[iwrk], &liwrk, &rwork[irwork], &ierr );        #elif TREVC_VERSION == 3        magma_ztrevc3( side, MagmaBacktransVec, select, n, A, lda, VL, ldvl,                       VR, ldvr, n, &nout, &work[iwrk], liwrk, &rwork[irwork], &ierr );        #elif TREVC_VERSION == 4        magma_ztrevc3_mt( side, MagmaBacktransVec, select, n, A, lda, VL, ldvl,                          VR, ldvr, n, &nout, &work[iwrk], liwrk, &rwork[irwork], &ierr );        #elif TREVC_VERSION == 5        magma_ztrevc3_mt_gpu( side, MagmaBacktransVec, select, n, A, lda, VL, ldvl,                              VR, ldvr, n, &nout, &work[iwrk], liwrk, &rwork[irwork], &ierr );        #else        #error Unknown TREVC_VERSION        #endif    }    time_sum += timer_stop( time_trevc );    flop_sum += flops_stop( flop_trevc );    if (wantvl) {        /* Undo balancing of left eigenvectors         * (CWorkspace: none)         * (RWorkspace: need N) */        lapackf77_zgebak( "B", "L", &n, &ilo, &ihi, &rwork[ibal], &n,                          VL, &ldvl, &ierr );        /* Normalize left eigenvectors and make largest component real */        for (i = 0; i < n; ++i) {            scl = 1. / magma_cblas_dznrm2( n, VL(0,i), 1 );            blasf77_zdscal( &n, &scl, VL(0,i), &ione );            for (k = 0; k < n; ++k) {                /* Computing 2nd power */                d__1 = MAGMA_Z_REAL( *VL(k,i) );                d__2 = MAGMA_Z_IMAG( *VL(k,i) );                rwork[irwork + k] = d__1*d__1 + d__2*d__2;            }            k = blasf77_idamax( &n, &rwork[irwork], &ione ) - 1;  // subtract 1; k is 0-based            tmp = MAGMA_Z_CNJG( *VL(k,i) ) / magma_dsqrt( rwork[irwork + k] );            blasf77_zscal( &n, &tmp, VL(0,i), &ione );            *VL(k,i) = MAGMA_Z_MAKE( MAGMA_Z_REAL( *VL(k,i) ), 0 );        }    }    if (wantvr) {        /* Undo balancing of right eigenvectors         * (CWorkspace: none)         * (RWorkspace: need N) */        lapackf77_zgebak( "B", "R", &n, &ilo, &ihi, &rwork[ibal], &n,                          VR, &ldvr, &ierr );        /* Normalize right eigenvectors and make largest component real */        for (i = 0; i < n; ++i) {            scl = 1. / magma_cblas_dznrm2( n, VR(0,i), 1 );            blasf77_zdscal( &n, &scl, VR(0,i), &ione );            for (k = 0; k < n; ++k) {                /* Computing 2nd power */                d__1 = MAGMA_Z_REAL( *VR(k,i) );                d__2 = MAGMA_Z_IMAG( *VR(k,i) );                rwork[irwork + k] = d__1*d__1 + d__2*d__2;            }            k = blasf77_idamax( &n, &rwork[irwork], &ione ) - 1;  // subtract 1; k is 0-based            tmp = MAGMA_Z_CNJG( *VR(k,i) ) / magma_dsqrt( rwork[irwork + k] );            blasf77_zscal( &n, &tmp, VR(0,i), &ione );            *VR(k,i) = MAGMA_Z_MAKE( MAGMA_Z_REAL( *VR(k,i) ), 0 );        }    }CLEANUP:    /* Undo scaling if necessary */    if (scalea) {        // converged eigenvalues, stored in WR[i+1:n] and WI[i+1:n] for i = INFO        magma_int_t nval = n - (*info);        magma_int_t ld   = max( nval, 1 );        lapackf77_zlascl( "G", &izero, &izero, &cscale, &anrm, &nval, &ione, w + (*info), &ld, &ierr );        if (*info > 0) {            // first ilo columns were already upper triangular,            // so the corresponding eigenvalues are also valid.            nval = ilo - 1;            lapackf77_zlascl( "G", &izero, &izero, &cscale, &anrm, &nval, &ione, w, &n, &ierr );        }    }    #if defined(VERSION3)    magma_free( dT );    #endif        timer_stop( time_total );    flops_stop( flop_total );    timer_printf( "dgeev times n %5d, gehrd %7.3f, unghr %7.3f, hseqr %7.3f, trevc %7.3f, total %7.3f, sum %7.3f/n",                  (int) n, time_gehrd, time_unghr, time_hseqr, time_trevc, time_total, time_sum );    timer_printf( "dgeev flops n %5d, gehrd %7lld, unghr %7lld, hseqr %7lld, trevc %7lld, total %7lld, sum %7lld/n",                  (int) n, flop_gehrd, flop_unghr, flop_hseqr, flop_trevc, flop_total, flop_sum );    work[0] = MAGMA_Z_MAKE( (double) optwrk, 0. );    return *info;} /* magma_zgeev */

//.........这里部分代码省略.........        #endif        lapackf77_cheevd(jobz_, uplo_,                         &n, A, &lda,                         w, work, &lwork,                         #ifdef COMPLEX                         rwork, &lrwork,                         #endif                         iwork, &liwork, info);        return *info;    }    /* Get machine constants. */    safmin = lapackf77_slamch("Safe minimum");    eps = lapackf77_slamch("Precision");    smlnum = safmin / eps;    bignum = 1. / smlnum;    rmin = magma_ssqrt(smlnum);    rmax = magma_ssqrt(bignum);    /* Scale matrix to allowable range, if necessary. */    anrm = lapackf77_clanhe("M", uplo_, &n, A, &lda, rwork);    iscale = 0;    if (anrm > 0. && anrm < rmin) {        iscale = 1;        sigma = rmin / anrm;    } else if (anrm > rmax) {        iscale = 1;        sigma = rmax / anrm;    }    if (iscale == 1) {        lapackf77_clascl(uplo_, &izero, &izero, &d_one, &sigma, &n, &n, A,                         &lda, info);    }    /* Call CHETRD to reduce Hermitian matrix to tridiagonal form. */    inde = 0;    indtau = 0;    indwrk = indtau + n;    indrwk = inde + n;    indwk2 = indwrk + n * n;    llwork = lwork - indwrk;    llwrk2 = lwork - indwk2;    llrwk = lrwork - indrwk;    magma_timer_t time=0;    timer_start( time );    magma_chetrd_mgpu(ngpu, 1, uplo, n, A, lda, w, &rwork[inde],                      &work[indtau], &work[indwrk], llwork, &iinfo);    timer_stop( time );    timer_printf( "time chetrd = %6.2f/n", time );    /* For eigenvalues only, call SSTERF.  For eigenvectors, first call       CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the       tridiagonal matrix, then call CUNMTR to multiply it to the Householder       transformations represented as Householder vectors in A. */    if (! wantz) {        lapackf77_ssterf(&n, w, &rwork[inde], info);        magma_smove_eig(range, n, w, &il, &iu, vl, vu, m);    }    else {        timer_start( time );        magma_cstedx_m(ngpu, range, n, vl, vu, il, iu, w, &rwork[inde],                       &work[indwrk], n, &rwork[indrwk],                       llrwk, iwork, liwork, info);        timer_stop( time );        timer_printf( "time cstedc = %6.2f/n", time );        timer_start( time );        magma_smove_eig(range, n, w, &il, &iu, vl, vu, m);        magma_cunmtr_m(ngpu, MagmaLeft, uplo, MagmaNoTrans, n, *m, A, lda, &work[indtau],                       &work[indwrk + n * (il-1)], n, &work[indwk2], llwrk2, &iinfo);        lapackf77_clacpy("A", &n, m, &work[indwrk + n * (il-1)], &n, A, &lda);                timer_stop( time );        timer_printf( "time cunmtr + copy = %6.2f/n", time );    }    /* If matrix was scaled, then rescale eigenvalues appropriately. */    if (iscale == 1) {        if (*info == 0) {            imax = n;        } else {            imax = *info - 1;        }        d__1 = 1. / sigma;        blasf77_sscal(&imax, &d__1, w, &ione);    }    work[0]  = magma_cmake_lwork( lwmin );    rwork[0] = magma_smake_lwork( lrwmin );    iwork[0] = liwmin;    return *info;} /* magma_cheevd_m */